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This  material  behavior  was  also  predicted  using  the  Taylor-type  rate  dependent  pxjlycrystal 
model  by  Asaro  and  Needleman.  Simulations  of  the  experiments  were  conducted  in 
personal  computers,  using  a  recently  developed,  highly  efficient,  fully  implicit,  time 
integration  scheme  by  Kalidindi  et  al.  In  the  initial  phase  of  the  simulation,  the  evolution  of 
ii'ic.  constituent  single  ciystal  slip  system  dciormation  resistance  was  estimated  using  the 
experimentally  determined  compressive  stress-strain  behavior  of  the  polycrystal.  With  this 
crystal  constitutive  behavior,  the  stress-strain  behavior  of  the  polycrystal  for  simple  shear 
was  computed.  In  addition,  the  evolution  of  the  crystallographic  texture  was  computed  for 
both  compression  and  simple  shear  tests.  Both  the  shear  and  axial  stress-strain  behaviors 
for  the  simple  shear  test  compared  reasonably  well  with  the  experimental  results. 
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INTRODUCTION 


In  an  event  such  as  the  penetration  of  an  armor,  the  armor  material  undergoes  a 
large  deformation  at  a  high  strain  rate  and  elevated  temperature  before  failure.  When 
numerically  simulating  such  a  rate  dependent  finite  deformation  event,  predictions 
significantly  depend  on  the  constitutive  models  used  to  represent  the  plastic  deformation 
behavior  of  the  material.  Most  computer  codes  contain  constitutive  models  that  are 
extensions  of  the  small  deformation  representation.  In  most  of  these  codes,  extension  to 
multiaxiality  is  achieved  through  the  assumption  that  there  exists  a  unique  effective  stress- 
strain  behavior  for  a  particular  material  at  a  given  strain  rate  and  temperature.  Finite  strain 
deformation  experiments  show  this  assumption  is  not  valid  for  large  deformations.  That  is, 
tension,  compression,  and  torsion  experimental  data  for  large  deformation  cannot  be 
correlated  by  an  effective  stress-strain  representation  ^  (also  from  the  experiments  repoitef'' 
in  this  report).  The  torsional  eticctive  stress-strain  curve  lies  significantly  below  both  the 
tensile  and  the  compressive  stress-strain  curves  at  large  strains.  Therefore  to  obtain  a  better 
prediction  of  the  deformation  at  finite  strain,  it  is  necessary  that  the  constitutive  models  in 
computer  codes  follow  the  response  of  the  material  at  finite  deformation  more  accurately. 

To  improve  upon  the  extension  of  the  small  deformation  classical  model  (Prandtl- 
Reuss),  finite-deformation  behavior  has  been  modeled  in  two  different  ways:  1) 
Macroscopic  phenomenological  modeling.  Here,  functional  forms  for  plastic  stretching, 
plastic  spin  and  back  stress  tensors  are  assumed.  2)  Microscopic  (polycry sta*'' 
phenomenological  modeling.  In  this  case,  the  phenomenological  representation  is  of  tf  , 
single  crystal  slip  behavior,  the  macroscopic  plastic  constitutive  behavior  is  derived  from  it. 
Also  the  evolution  of  the  crystallographic  texture  (crystallographic  anisotropy)  can  be 
deduced  from  this  type  of  formulation.  In  addition,  with  the  insight  into  the  micro¬ 
mechanics  of  the  plastic  deformation  process,  the  polycrystal  modeling  can  lead  us  to  select 
functional  forms  for  macroscopic  modeling  that  will  better  represent  the  physics  of 
microscopic  plastic  deformation. 

In  the  polycrystal  phenomenological  plasticity  modeling,  three  methods  have  been 
used  to  deduce  the  polycrystal  constitutive  behavior  from  the  deformation  of  the  individual 
crystals  (grains):  1)  The  Taylor  method,^  where  each  grain  is  assumed  to  deform  with  the 
same  deformation  gradient  as  the  polycrystal,  thus  satisfying  compatibility  among  the 
grains,  but  violating  the  stress  equilibrium  across  the  grains.  2)  Sachs  method,^  where 
each  grain  is  assumed  to  be  subjected  to  the  same  stress  tensor  as  the  polycrystal,  thus 
satisfying  the  stress  equilibrium  condition,  but  violating  the  compatibility  across  the  grains. 
3)  Self-consistent  method,  where  each  grain  is  considered  as  an  inclusion  embedded  in  a 
homogeneous  infinitely  extended  matrix.^-^A  Unfortunately,  there  are  not  that  many 
experimental  evaluation  of  these  polycrystal  schemes  in  the  literature. 

1.  SEMIATIN,  S.  L.,  LAHOTl,  G.  D.  AND  JONAS,  J.  J.  Mechanical  Testing,  Metals  Handbook, 
Ninth  Editon,  v.  8,  p.  164. 

2.  TAYLOR,  G.  I.  J.  Plastic  Strain  in  Metals.  In.st.  Mctal.s,  v.  62,  1938,  p.  307. 

3.  SACHS,  G.,  Z.  Vcr.  Dtsch.  Ing.,  v.  72,  1928,  p.  734. 

4.  KRONER,  E.  Berechning  der  elastischen  Konstanten  des  Vielkristalls  aus  den  Konstanten  des 
Pinkristalls.T,.  Phys.,  v.  151,  1958,  p.504-518. 

5.  nUDIANSKY,  B.  AND  WU,  T.  T.  Theoretical  prediction  of  plastic  strains  of  poly  cry. stals.  In  Proc. 
of  the  4th  U.  S.  National  Congrc-ss  of  Applied  Mechanics,  ASME,  New  York,  1962,  p.  1 175-1185. 

6.  1111, L,  R.  Continuum  micro-mechanics  of  elastopla.stic  polycrystals.  J.  Mech.  Phys.  Solids,  v.  13, 
1965,  p  89-101. 
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For  elastic-viscoplastic  finite  deformation  via  slip  by  dislocation  glide,  Asaro  and 
Needleinan^  have  proposed  a  Taylor-typc  model  to  represent  the  polycrystai  behavior  from 
its  individual  crystals.  Harren  ct  al.®  have  conducted  a  detailed  analysts  of  the  finite 
deformation  under  shear  (torsion)  using  this  model.  They  have  observed  that  the  axial 
behavior  of  a  single  phase  polycrystalline  material  in  a  torsion  test  is  very  sensitive  to  strain 
hardening,  strain  rate  sensitivity,  and  latent  hardening  behavior  of  the  constituent  single 
crystal.  Therefore,  finite  deformation  torsion  testing  can  be  used  to  evaluate  and  calibrate 
polycrystal  plastic  models  as  well  as  macroscopic  m^els. 

Four  types  of  experiments  are  commonly  used  to  evaluate  plastic  constitutive 
behavior  of  materials:  1)  uniaxial  tensile  testing,  2)  uniaxial  compression  testing,  3)  thin 
walled  torsion  testing,  and  4)  multiaxial  testing.  With  tensile  testing,  it  is  not  possible  to 
obtain  experimental  data  for  large  deformations  due  to  the  necking  instability. 
Compression  testing  can  be  used  to  obtain  large  deformation  data  up  to  strains  of 
approximately  1.5.  Though  multiaxial  testing  is  useful  to  study  constitutive  models,  it  is 
complex  to  perform,  especially  to  large  strains.  However,  with  torsion  testing,  it  is 
possible  to  obtain  experimental  finite  deformation  data  up  to  very  large  strains.  In  addition, 
axial  stress  or  strain  measurement  from  the  torsion  tests,  depening  on  whether  these  tests 
are  conducted  under  fixed-end  or  free-end  conditions,  is  a  very  sensitive  measure  and  can 
be  used  to  evaluate  both  polycrystal  and  macroscopic  finite  deformation  constitutive 
models.  However,  there  is  a  lack  of  torsion  test  data  with  axial  stress  or  strain 
measurements  in  the  literature. 

A  set  of  finite  deformation  experiments  (fixed  and  free-end  thin  wall  torsion  and 
uniaxial  compression)  were  conducted  for  the  experimental  evaluation  of  the  finite 
deformation  models.  This  report  gives  an  evaluation  of  the  Taylor-type  rate  dependent 
model  proposed  by  Asaro  and  Needleman^  using  the  experimental  stress-strain  behavior 
from  these  experiments.  The  predicted  textures  from  the  rate  dependent  Taylor  model  are 
also  given  in  this  report.  Comparison  of  the  predicted  textures  with  the  experimental 
textures  will  be  reported  in  a  future  publication. 


THE  TAYLOR-TYPE  POLYCRYSTAL  MODEL 

An  overview  of  the  Taylor-type  rate  dependent  polycrystai  model  by  Asaro  and 
Needleman  is  given  in  this  section.  In  our  work,  we  followed  the  formulation  of  the  time 
integration  of  the  polycrystai  model  for  various  boundary  conditions  as  proposed  by 
Kalidindi  et  al.^  The  summary  of  their  integration  algorithm  is  also  given  in  this  section. 


7.  ASARO,  R.  J.  AND  NEEDLEMAN,  A.  Texture  Development  and  Strain  Hardening  in  Kate 
Dependent  I’olye.rystah.  Acta  Mclall.,  v.  33, 1985,  p.  923. 

8.  HARREN,  S.,  LOWE,  T.  C,  ASARO,  R.  J.  AND  NEEDLEMAN,  A.  Analysis  of  large-strain  shear 
in  rate-dependent  face-centred  cubic  polycrystals:  correlation  of  micro-  and  macromechanics,  Phil. 
Trans.  R.  Soc.,  v.  328,  1989,  p.  443-500. 

9.  KALIDINDI,  S.  R.,  BRONKHORST,  C.  A.  AND  ANAND,  L.  Crystallographic  Texture  Evolution 
in  Bulk  Deformation  Processing  of  FCC  Metals,  to  appear  in  Journal  of  the  Mechanics  and  Physics  of 
Solids. 


With  this  formulation,  the  problem  could  be  solved  with  PCs  due  to  the  efficiency  of  the 
algorithm.  The  disadvantage  of  this  formulation  is  that  it  can  be  applied  only  with  the 
prescribed  displacement  boundary  conditions. 

Kinematics  of  the  Single  Crystal 

Let  F  be  the  global  deformation  gradient  of  the  polycrystalline  material.  With  the 
Taylor  assumption,  then  the  deformation  gradient  of  each  crystal  of  the  polycrystal  will 
have  the  same  deformation  gradient  F.  By  multiplicative  decomposition  of  the  deformation 
gradient  into  elastic  and  plastic  components,  we  obtain: 

F  =  FeFP  (1) 

where  F®  is  the  elastic  and  FP  is  the  plastic  parts  of  the  deformation  gradient.  Here  FP  is 
the  deformation  gradient  due  to  plastic  slip  only  fdet  FP=  1);  it  represents  an  elastically 
unloaded  intermediate  relaxed  configuration  of  a  stress  free  state  with  the  orientation  of  the 
crystal  lattice  coinciding  with  the  crystal  lattice  in  the  reference  configuration.  The 
deformation  gradient  F®  represents  both  rotation  and  elastic  stretching  of  the  crystal  lattice 
(det  F®>1).  Then  the  velocity  gradient  L: 


L  =  ^  F-1  =  f"®  F®-1  +  F®  LB  F®-1 

(2) 

where 

LB  =  fP  FP-l 

(3a) 

or 

^P  =  LB  FP 

(3b) 

Here  the  notation  (°)  indicates  the  derivative  with  respect  to  time.  Since  LB  is  the  plastic 
velocity  gradient  in  the  relaxed  configuration,  LB  is  given  by 

12 

L8=I^S?j,  S«=sB®nB  (3c) 

a=l 


where  sB  and  ng  are  orthonormal  unit  vectors  representing  slip  directions  and  slip  planes, 
respectively,  for  each  slip  system  a  in  the  reference  configuration  and  the  is  the  plastic 

shearing  rate  in  the  ath  slip  system.  For  an  fee  crystal  there  are  twelve  slip  systems 
derived  from  its  four  (111)  slip  planes  and  three  <1 10>  slip  directions  for  each  slip  plane. 
Consequently,  unit  vectors  nB  represent  four  (111)  planes  and  the  unit  vectors  sB 

represent  the  corresponding  three  <1 10>  directions  of  each  ath  slip  system  (in  the  Table  of 
section  Model  Predictions).  From  Equation  (2), 

12 

LP  =  DP  +  WP  =  F®  LB  F®-1  =  I  ^  S«,  S«  =  (F®  sB  )  <S>  (nB  F®-!)  (4a) 

a=l 

where  DP  and  WP  are  plastic  stretch  and  spin  tensors,  respectively,  and  are  given  by 
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12  12 

DP=  I  ^  sytn{S“),  and  WP=  I  asym{S«)  (4b) 
a=l  0=1 

The  rate  of  plastic  stress  power  per  unit  reference  volume  is  (see  Anand'®) 

&P  =  (Ce  T*):LP,  with  =  peT  pe  (5) 

and  T*  =  F*- 1  { (det  Pe)  T }  pe-T.  (6) 

The  stress  T  is  Cauchy  stress  of  the  crystal.  The  second  Piola-Kirchhoff  stress,  T*,  is  the 
elastic  work  conjugate  of  the  Lagrangean  elastic  strain  measure 

pe  =  (1/2)  {Ce  -  1),  with  C®  =  peT  pc  (7) 

But  the  plastic  stress  power  per  unit  volume  is 

12 

^  (8) 

0=1 

where  x“  is  the  resolved  shear  stress  in  the  oth  slip  system.  Therefore,  from  Equations  (5) 
and  (8)  resolved  shear  stress 

-ja  =  (ce  T*):Sg  (9) 

Crystal  Plasticitv  Constitutive  Model 

The  relationship  between  the  resolved  shear  stress,  and  resolved  shear  strain 
rate,  Yo  ,  is  given  by  the  viscoplastic  rate  power  law 

^  =  Yo  I  ^“/s”  sign(T“)  ( 1 0) 

where  Yo  is  the  reference  shear  strain  rate,  s“  is  the  resistance  of  the  slip  system  o  and  m  is 
the  strain  rate  sensitivity  parameter  (assumed  to  be  the  same  for  all  slip  systems).  Also,  the 

slip  resistance  s“  can  be  considered  as  the  slip  resistance  of  the  slip  system  o  at  the 
reference  shearing  rate,  % 

The  slip  resistance  s“  is  evolved  according  to 


10.  ANAND,  L.  Constitutive  Equations  for  Hot-Working  of  Metals.  Int  J.  Plast.,  v.  1,  1985,  p.  213- 
231. 
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(11) 


12 

s“=  1^1 

p=l 


where  h“P  is  the  instantaneous  hardening  modulus  of  the  slip  system  a  due  to  a  shearing 

on  the  slip  system  p.  The  form  of  the  hardening  moduli  (Asaro  and  Needleman)  are  given 
by 


h“P  =  h(sP)  q“P  (no  sum  of  P)  (12) 

where  q“P  is  a  12  by  12  matrix  describing  the  relationship  between  self  and  latent 
hardening  behavior  of  the  crystal  slip  systems.  The  ratio  of  latent  hardening  to  self 
hardening  is  assumed  to  be  unity  for  coplanar  and  a  value  of  q  for  non-coplanar  systems. 
From  Kalidindi  et  al.  the  single  slip  hardening  modulus,  h(s),  is  assumed  to  be  given  by 

h(s)  =  ho(  l-f)a  (13) 

*s 

where  constants  ho,  Sg  and  a  are  assumed  to  be  the  same  for  all  slip  systems.  Here,  Sg 
represents  the  saturated  slip  resistance  and  ho  is  a  measure  of  the  initial  hardening  modulus 
of  the  single  slip  system.  From  this  assumed  single  slip  hardening  behavior 

for  a  #  -1,  where  So  is  the  initial  slip  resistance  and  the  I  y  I  is  the  absolute  value  of  shear. 
Here,  the  Bauschinger  effect  is  assumed  to  be  negligible. 

From  finite  elasticity,  the  reduced  constitutive  equation  for  an  elastic  material  in  a 
grain  can  be  written  in  the  form 

T*  =  A'E‘5]  (I'ii 

where  L  is  the  fourth  order  elastic  tensor.  Assuming  all  the  grams  are  of  equal  volume, 
with  the  Taylor  assumption  that  the  local  deformation  gradient  of  each  grain  is  the  same  as 
the  global  deformation  gradient,  then  the  volume  average  Cauchy  stress  T^  is 

N 

(16) 

^  k=l 


where  T*^  is  the  Cauchy  stress  of  each  grain  and  N  is  the  number  of  grains.  The  volume 
average  stress  is  the  global  stress  response  of  the  polycrystal  (Asaro  and  Needleman). 
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In  the  reformulation  of  the  Asaro-Needleman  polycrystal  model  for  a  more  efficient 
tiMe-integration  scheme,  Kalidindi  et  al.^  have  assum^  the  following  approximations; 

For  small  elastic  stretches: 


5 

1 

+ 

<N 

III 

(17) 

and 

xa  =T*:S“ 

(18) 

Numerical  Procedure  for  Model  Predictions 

In  the  simulation  of  the  experiments  reported  in  this  report,  the  time-integration 
procedure  proposed  by  Kalidindi  et  al.  was  used  to  integrate  the  above  polycrystal  model. 
This  formulation  of  the  time-integration  procedure  is  summarized  in  this  section. 

The  problem  of  integrating  the  above  polycrystal  model  can  be  summarized  as 
follows.  Let  T  be  the  time  At  later  than  the  current  time  t.  Then  the  problem  is  to  find  the 
list  {FP(t),  s®(t),  T(t)}  for  each  grain  at  time  x  when  the  current  state  {FP(t),s“(t),  T(t)} 
of  all  the  grains  and  the  global  deformation  gradients  at  times  t  and  t  are  given.  Once  the 
stress  at  each  grain  is  known,  global  stress  behavior  at  x  can  be  calculated  from  Equation 
(16).  The  new  texture  at  time  x  for  each  crystal  can  be  computed  from  (see  Equation  4a) 


s'?  =  Fc(x)  s8f 

(19) 

ni?  =  Fe-T(x)  ng 

(20) 

where  (s'?,  n'?)  is  the  slip  system  in  the  deformed  configuration.  Since  the  slip  systems  in 
the  undeformed  configuration  (sgf,  n^)  are  orthonormal,  the  slip  systems  in  the  deformed 
configuration  are  also  orthonormal.  If  F®(x)  is  found,  then  the  slip  system  in  the  deformed 

f'onfiguration  for  each  crystal  at  the  time  x  can  be  computed  and  hence  the  pole  figure  can 
be  constructed  for  the  polycrystal. 

The  starting  point  for  the  formulation  of  the  time  integration  procedure  is  the  result 
obtained  by  Weber  and  Anand’*  in  the  implicit  time-integration  of  the  evolution  equation  for 
f'P  given  in  the  Equation  (3b): 


FP(x)  =  exp  {At  LP(x)}  FP(t) 


(21) 


1 1.  WEBER,  G.  AND  ANAND,  L.  Finite  Deformation  Conslilulivc  Equations  and  a  Time  Integration 
Procedure  for  Isotropic,  Hypcrelastic-Viscoplaslic  Solids.  Comp.  Methods  Appl.  Mcch.  Eng.,  v.  79, 
1990,  p.  173-202. 
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By  neglecting  higher  order  terms  of  Ay®  for  small  A*j®  and  using  Equation  (3c) 


or 


12 

FP(t)=  {1+  FP(t) 

a=i 


12 

FP-l(t)  =  FP->(t)  {I- 

a=l 


where  Ay®  -  "^(x)  At. 

From  Equations  (23),  (18)  and  (10)  A7®  is  given  by 


(22a) 

(22b) 

(23) 


Ay®  =  Yo  sign(T*:S“)  (24) 

Substituting  Equations  (7)  and  (1)  in  the  elastic  constitutive  Equation  (15),  following 
equations  are  obtained  (see  Kalidindi  et  al."^  for  details): 


12 

T*(t)  =  T*®-  XaY“C®  (25) 

a=l 


12 

s®(x)  =  s®(t)  +  I  AyP  I 

(J=l 

where 


(26) 


T*®  =  ^A+^^(trA)  1-^1, 
A  =  FP-T(t)  FT(t)  F(x)  FP-l(t), 

C®  =  )t  B®  +  (trB®)  1, 

B®  =  A  Sg+SgTA, 


(27) 

(28) 

(29) 


Equations  (25)  and  (26)  are  two  sets  of  nonlinear  simultaneous  equations  (T*(t)  -  six 
unknowns  -  and  s®(x)  -  twelve  unknowns)  and  can  be  solved  using  a  modified  Newton- 
Raphson  type  algorithm  (see  Kalidindi  et  al."^  for  the  details  of  the  algorithm).  Once  T*(t) 
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and  s“(x)  are  known,  FP(t)  is  obtained  using  Equation  (22a)  and  F®(x)  is  obtained  from 
Equation  (1).  The  new  texture  at  time  x  can  be  found  by  the  polar  decomposition  of  F®(x) 

into  its  elastic  lattice  rotation  component  (R®(x))  and  elastic  stretching  component  and  by 
using  the  Equations  (19)  and  (20)  to  obtain: 


sSf  =  Re(x)  sg 
(31) 

n(?  =Re-T(x)ng.  (32) 

where,  s<f  (slip  direction)  and  (slip  plane)  are  unit  vectors  representing  the  slip  system 
a  at  time  x. 

A  computer  program  that  solves  this  nonlinear  problem  is  given  in  the  Appendix. 


EXPERIMENTAL  STRESS-STRAIN  BEHAVIOR 


Material 

Test  specimens  were  machined  from  oxygen-free-high-conductivity  (OFHC) 
copper  (99.99%  Copper).  OFHC  copper  was  receiv^  as  bar  stocks  in  the  work  hardened 
condition.  Before  machining,  the  copper  rods  were  annealed  at  400°C  in  argon  for  one 

hour.  This  annealing  treatment  produced  an  average  grain  diameter  of  45  )im. 

Experimental  Procedure 

Compression  Tests 

All  the  compression  tests  were  conducted  with  cylindrical  specimens  0.5  in  (12.7 
mm)  in  diameter  with  a  nominal  gage  length  of  0.5  in  (12.7  mm)  using  an  Instron  screw 
driven  testing  machine.  The  specimen’s  ends  were  grooved  as  shown  in  the  Figure  1  to 
retain  the  lubricant  for  a  'onger  test  time.  To  record  the  initial  5%  strain,  all  the  specimens 
were  strain  gaged.  The  ends  of  the  specimens  were  lubricated  with  thin  teflon  tapes.  The 
relative  displacement  of  the  ends  of  the  specimens  were  monitored  using  a  jig  \yith  a  DCDT 
transducer.  All  the  tests  were  conducted  at  a  nominal  engineering  strain  rate  of  5.7  x  lO'^ 
s'l.  The  stress  strain  curve  for  engineering  strains  greater  than  5%  was  obtained  from  the 
DCDT  information  and  using  the  assumpdon  that  the  deformation  is  incompressible.  The 
discontinuity  observed  at  the  5%  strain  level  between  the  two  curves  was  reconciled  by 
moving  the  starting  point  of  the  stress-strain  curve  for  strains  greater  than  5%  to  coincide 
with  the  end  point  of  the  stress-strain  curve  for  strains  less  than  5%.  This  discrepancy 
arises  due  to  the  presence  of  the  grooves  at  the  loading  interfaces  of  the  specimen. 
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Compression  Test  Results 


Stress-strain  results  for  the  compression  tests  are  given  in  Figure  2.  These  data 
were  obtained  using  five  different  specimens.  The  maximum  strain  that  was  obtained  was 


1.30. 


L  -  0.500- 
0  ■  0.500' 


Figure  1.  Compression  Specimen 
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-I>  TEST-CS 
-B-  TEST.C6 
-B*  TEST.C3 
-o  TEST-C14 
-O  TEST-C9 
— SIMULATION 


0.6  0.8 

TRUE  STRAIN 


Figure  2.  Experimental  Uniaxial  Compression  Stress-Strain  Data. 
Simulated  Compression  Stress-Strain  Curve  (for  the  Single  Slip  Behavior 
So  =  16  MPa;  ho  =  700  MPa;  Ss  =  1.55  MPa;  a  =  3.8)  is  also  given. 


Torsion  Tests 


Specimen  geometry  of  the  torsion  test  specimen  is  given  in  Figure  3.  This 
specimen  is  a  modified  Lindholm-ty|3e  specimen  that  was  used  by  White  et  alJ^  The  gage 
portion  of  the  torsion  specimen  is  a  thin  wall  tube  of  0.240  in  (6.096mm)  length  and 
external  and  internal  diameters  0.810  in  (20.574mm)  and  0.750  in  (19.050mm) 
respectively  (0.030  in  wall  thickness).  The  torsion  tests  were  conducted  using  both  free 
and  fixed-end  conditions  at  a  shear  strain  rate  of  0.001  s'^  For  the  tests  that  were 
conducted  under  fixed-end  conditions,  the  axial  strain  was  assumed  to  be  zero.  Though  the 
axial  strain  in  the  gage  sections  were  not  directly  measured,  axial  extension  just  outside  the 
gage  section  was  monitored  for  one  free-end  and  one  fixed-end  test  using  a  specially 
designed  axial  extensometer  (see  Figure  3).  All  the  torsion  tests  were  conducted  using  a 
tension-torsion  hydraulic  test  machine.  Specimens  were  attached  to  the  test  machine  using 
a  pair  of  hydraulic  collet  grips. 

Shear  stress  was  derived  using  the  expression 

Shear  stress  =  Torque  /  (27crm^t)  (33) 

where  rm  is  the  mean  radius  of  the  gage  section  and  t  is  the  wall  thickness  of  the  gage  area. 
The  axid  stress  was  obtained  by 

Axial  stress  =  Axial  Force  /  (7i(R2  -  r^))  (34) 

where  R  is  the  outer  radius  and  r  is  the  inner  radius.  Axial  strain  is  given  by  the  expression 

Axial  Strain  =  Axial  Extension/Gage  Length  (35) 

and  the  shear  strain  is  given  by  the  expression 

Shear  Strain  =  (R0)/Gage  Length  (36) 

where  0  (in  radians)  is  the  relative  rotation  at  the  grips. 

Torsion  Test  Results 

Stress-engineering  shear  strain  curves  obtained  for  fixed-end  torsion  tests  are  given 
in  Figure  4.  The  axial  stress  is  initially  compressive  and  goes  through  a  minimum  at  an 
approximate  shear  strain  of  1  with  the  monotonic  increase  of  shear  strain;  as  the  shear  strain 
is  increased  further,  the  axial  stress  becomes  tensile  at  an  approximate  strain  of  2.25.  The 
axial  strain  that  was  obtained  by  monitoring  the  displacement  just  outside  the  end  of  the 
shoulder  section  is  also  given  in  the  Figure  4.  The  axial  strain  monitored  this  way  was 


12.  WHITE,  C.  S.,  BRONKHORST,  C.  A.  AND  ANAND,  L.  An  Improved  fsotropic-Kinemaiic 
Hardening  Model  for  Moderate  Deformation  Metal  Plasticity.  Mech.  MaL,  v.  10,  1990,  p.  127-147. 
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approximately  zero  up  to  a  shear  strain  of  1.75  and  increases  after  this  strain  (positive  strain 
indicates  an  extension  in  the  gage  section). 

Shear  stress-shear  strain  behavior  obtained  using  ffee-end  tests  is  given  in  Figure  5. 
Also  in  this  figure,  the  axial  strain  behavior  is  given  as  a  function  of  the  shear  strain.  Axial 
strains  were  obtained  both  from  the  machine  and  from  the  extensometer.  As  expected  the 
axial  strain  measured  from  the  extensometer  was  lower  than  measured  by  the  machine. 
Here,  positive  strain  indicates  an  extension.  The  axial  strain  increases  and  remains 
approximately  constant  after  a  shear  strain  of  about  1. 

Both  fixed-end  and  free-end  stress-strain  behavior  is  given  in  Figure  6.  The  shear 
stress-shear  strain  behavior  is  approximately  equal  and  therefore  insensitive  to  the  method 
of  testing  -  i.e.,  whether  the  type  of  testing  is  fixed-end  or  free-end. 


i 


Figure  3.  Torsion  Specimen  with  the  Extensometer 
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STRESS  (MPa)  /  STRAIN  (per  1000) 


3 


SHEAR  STRESS  (MPa)  /  STRAIN  (per  1000) 


Figure  5.  Shear  Stress/Axial  Strain  versus  Shear  Strain  From  Free-End 

Torsion  Tests 
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SHEAR  STRESS  (MPa) 


Figure  6.  Comparison  of  Experimental  Fixed-End  and  Free-End  Torsion 

Data 
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MODEL  PREDICTIONS 


The  computer  program  given  in  the  Appendix  was  used  to  obtain  the  polycrystal 
model  predictions  of  the  uniaxial  compression  and  fixed-end  torsion  (simple  shear)  tests 
and  the  evolution  of  crystallographic  texture.  Following  Kalidindi  et  al.,  the  strain  rate 
sensitivity  parameter,  m,  for  OFHC  copper  was  taken  as  0.012.  Latent  hardening,  q,  for 
the  copper  was  assumed  to  be  1.4  for  the  computation.  Also  simulations  were  conducted 
with  m  =  0.04;  q  =  1.4,  and  m  =  0.012;  q  =  1.0  to  evaluate  the  effect  of  higher  strain  rate 
sensitivity  and  lower  latent  hardening. 

Computations  were  conducted  using  a  16  MHz  IBM  compatible  80386  computer 
with  80387  floating  point  coprocessor  support  or  a  Mac  Ilex  computer.  The  initial 
crystallographic  texture  was  assumed  to  be  isotropic  (the  grains  were  randomly  oriented); 
the  initial  distribution  of  the  crystals  were  taken  from  Molinari  et  al.'^  and  was  given  as  300 
sets  of  Euler  angles.  Each  set  of  Euler  angles  would  rotate  the  axes  of  the  corresponding 
crystal  to  coincide  with  the  global  fixed  rectangular  Cartesian  axes.  Schmid  tensor,  Sgf, 
was  computed  with  respect  to  the  global  coordinate  system  for  each  crystal  by 

S?i=QS?QT  (37) 

where  S?  is  defined  with  respect  to  an  orthonormal  basis  associated  with  the  crystal  lattice 
<100>  directions.  Twelve  slip  systems  which  defines  S?  for  fee  crystal  are  given  in  the 
Table  1. 

The  onhogonal  matrix  Q  which  rotates  the  crystal  basis  to  coincide  with  the  global 
fixed  basis  is  given  by 


Q  = 


cos  <t>  sin  6  sin  <|>  cos  co 

-sin  (j)  sin  o)  cos  9  +cos  <t)  sin  to  cos  0 

-cos  <|)  sin  to  -sin  (j)  sin  to 

-sin  6  sin  to  cos  0  cos  6  cos  to  cos  9 


sin  ({)  sin  0 


-cos  <j>  sin  0 


sin  to  sin  0 


cos  to  sin  0 


cos  0 


(38) 


where  {0  <  6  <  27c,  0  <  0  <  rt,  0  <  to  <  27c)  are  the  three  Euler  angles  that  represent  this 
transformation.  Most  of  the  simulations  were  conducted  with  a  subset  of  100  crystals  after 
it  was  determined  that  there  was  no  appreciable  difference  in  the  results  whether  300  or  100 
c^stals  were  used  (see  sub-section  Fixed-end  Torsion  Behavior).  Both  of  these  crystal 
distributions  are  given  in  Figure  7  as  equal  area  projection  (111)  pole  figures. 


13.  MOLINARI,  A.,  CANOVA,  G.  R.  AND  AHZl,  S.  A  Seif-Consistent  Approach  of  the  Large 
Deformation  Polycrystal  Viscoplasticity.  Acta  Mctail.,  v.  35,  1987,  p.  2983-2994. 
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Figure  7.  {Ill}  Equal- Area  Projection  Pole  Figures  of  Initial  Randomly 

Oriented 

a).  ^00  and,  b).  300  Crystals 
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Table  1.  Twelve  slip  systems  for  fee  erystal 


slip  system 

a 

slip  plane 

n« 

slip  direction 

s“ 

1 

(  1  1  1)W3 

f-1  1  OI/V2 

2 

( 1  1  1)W3 

[-1 0 1W2 

3 

(  1  1  1)W3 

[  0-1  \]Hl 

4 

[  1  1  01/V2 

5 

(-1  1  1)/V3 

f  1 0 1W2 

6 

(-1  1  1)/V3 

[0-1 1W2 

7 

(  1-1  1)/V3 

[  1  1  0W2 

8 

(  1-1  1)W3 

[  0  1  11/V2 

9 

(  1-1  1)W3 

[  1  0-l]W2 

10 

( 1  1-1)/V3 

f  0 1 1W2 

11 

(  1  1-1)/a/3 

[  1 0 1W2 

12 

( 1  1-1)/V3 

[1-1  0]W2 

Uniaxial  Compression  Behavior 

The  time  integration  algorithm  diseussed  earlier  requires  that  the  deformation 
gradient  to  be  preseribed  for  the  deforming  body.  In  a  uniaxial  eompression  test,  axial 
strain  rate  is  known;  but  the  lateral  strain  rates  are  not  known.  Following  Bronkhorst  et 
al.,*^  for  the  uniaxial  eompression  test,  the  material  elements  are  assumed  to  follow  an 
isoehorie  motion  given  by 

X  =  exp{£t/2)Xi  ej  +  exp{^t/2)X2  ©2  exp{-^t)X3  63  (39) 


where  I  is  the  strain  rate,  X  {Xi,X2,X3}  and  x  are  the  original  and  eurrent  positions  of  a 
material  element  respeetively  and  ej  {i  =  1, 2,  3)  are  orthonormal  base  veetors  of  a  global 
rectangular  Cartesian  coordinate  system.  The  response  from  the  Taylor  model  for  this 
motion  is  the  volume  averaged  stress  deviator  T'  =  T  +  pi,  where  p  an  undetermined 
pressure.  In  a  uniaxial  compression  test,  because  the  lateral  tractions  are  zero,  the 
approximate  value  of  the  axial  Cauchy  stress  can  be  represented  by 

"^33  =  T'33  -  P-  where  p  =  (T 1 1  +  r22)/2  (40) 


14.  BRONKHORST,  C.  A..  KALIDINDl  S.  R.  AND  ANAND,  L.  Polycrystalline  Plasticity  and  the 
Evolution  of  Crystallographic  Texture  in  FCC  Metals,  to  be  published,  1991. 
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The  assumed  single  slip  hardening  modulus  and  resistance  of  the  constitutive 
behavior  of  the  single  crystal  is  given  by  the  Equations  (13)  and  (14).  Here,  the  unknown 
parameters  are  h^  -  the  initial  h^ness,  Sg  -  saturated  vdue  of  the  slip  resistance,  Sq  -  the 
initial  slip  resistance  or  the  initial  critical  resolved  shear  stress  (crss)  and  a  -  the  exponent. 
By  systematically  iterating  these  parameters,  the  calculated  macroscopic  stress-strain 
behavior  can  be  adjusted  to  represent  the  experimentally  observed  stress-strain  behavior. 
The  saturated  stress-strain  behavior  is  dependent  on  the  Sg  value,  initial  slope  of  the  stress- 
strain  curve  is  dependent  on  the  h^  value,  yield  point  is  dependent  on  the  Sq  value,  and  the 
exponent,  a,  determines  the  shape  of  the  stress-strain  curve  between  the  yield  point  and  the 
saturation  point.  By  systematically  changing  the  values  of  these  parameters,  a  simulated 
compression  stress-strain  curve  was  obtained,  which  approximately  represents  the 
experimental  stress-strain  behavior.  This  simulated  stress-strain  behavior  is  given  in 
Figure  2  with  the  experimental  stress-strain  curves.  Corresponding  parameters  of  the 
single  slip  behavior  are 

Sq  =  16  MPa,  ho  =  700  MPa,  Sj  =  155  MPa  and  a  =  3.8.  (41) 

These  parameters  for  single  slip  behavior  are  used  for  the  simulation  of  the  torsion  (simple 
shear)  stress-strain  behavior  in  the  next  section.  Thus,  by  curve  fitting  the  Taylor  model  to 
the  experimental  results  of  the  compression  tests,  the  parameters  for  simulating  the  torsion 
tests  are  fixed. 

For  this  single  slip  behavior,  simulated  stress-strain  behaviors  for  m  =  0.04,  q  = 
1.4  and  m  =  0.012,  q  =  1.0  are  given  in  the  Figure  8  for  comparisons.  As  the  strain  rate 
sensitivity  increases  or  latent  hardening  decreases,  for  uniaxial  compression  stress-strain 
loading,  the  simulation  predicts  a  lower  stress  for  the  same  strain. 

Fixed<encl  Torsion  (Simple  Shear)  Behavior 

Using  the  parameters  that  were  determined  from  the  simulation  of  the  compression 
tests,  fixed-end  torsion  tests  were  simulated.  In  the  fixed-end  torsion  tests,  the  material 
elements  are  assumed  to  follow  a  motion  described  by 

x  =  (Xi +7tX2)ei +X2e2  +  X3e3.  (42) 

where  y  is  the  engineering  shear  strain  rate.  The  simulated  shear  (T12),  axial  (T22)  and 
hoop  (Til)  stresses  as  functions  of  shear  strain  are  given  in  Figure  9  for  both  100  and  300 
crystals.  There  is  no  appreciable  difference  in  the  stress-shear  strain  behavior  for  100  or 
300  crystals  except  the  axial  and  hoop  stress  behavior  become  smoother  for  300  crystals. 
Therefore,  all  the  other  simulations  that  are  reported  in  this  report  were  carried  out  with  100 
crystals. 

Figure  10  gives  the  simulated  stress-shear  strain  behavior  (m  =  0.012  and  q  =  1.4) 
for  simple  shear  losing  with  the  corresponding  fixed-end  experimental  data  from  Figure  4. 
Simulated  axial  stresses  are  compressive  in  the  range  of  shear  strains  up  to  3  and  go 
through  a  compressive  maximum  stress  as  the  shear  strain  is  increased.  In  contrast,  the 
hoop  stress  is  tensile  and  increases  as  the  shear  strain  is  increased.  For  comparison  with 
this  simulated  stress-strain  behavior.  Figure  1 1  gives  the  calculated  stress-shear  strain 
behavior  for  a  higher  rate  sensitive  material  (m  =  0.04;  q  =  1.4)  and  for  a  lower  latent 
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STRESS  (MPa) 


hardening  (q  =  1.0;  m  =  0.012)  material.  In  both  of  these  simulations,  single  slip 
hardening  was  assumed  to  follow  the  behavior  given  by  the  parameters  in  (4 1 ). 


TRUE  STRAIN 


Figure  8.  Effect  of  Strain  Rate  Sensitivity  and  Latent  Hardening  on 
Simulated  Compression  Stress-Strain 
(So  =  16  MPa;  ho  =  700  MPa;  s*  =  155  MPa;  a  =  3.8) 
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Figure  9.  Comparison  of  Simulations  of  Simple  Shear  with  100  versus  300 

Crystals 
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STRESS  (MPa) 


SHEAR  STRAIN 


Figure  10.  Comparison  of  Fixed-End  Torsion  Experimental  and  Simulated 

Stress-Strain 
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STRESS  (MPa) 


200 


Figure  11.  Simulated  Simple  Shear  Stress-Strain  as  a  function  of 
Strain  Rate  Sensitivity  (m=  0.012  and  0.04)  and 
Latent  Hardening  (q  =  1.4  and  1.0) 

(single  slip  behavior  from  experimental  compression  data  with  m=0.012 

and  q=  1.4) 
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A  simulation  was  also  conducted  starting  from  compression  data  for  these  two 
cases:  1)  m  =  0.04;  q  =  1.4  -  a  higher  rate  sensitive  material,  for  which  single  slip  behavior 
was  found  to  be  given  by  the  Equation  41,  except  for  Ss  which  was  166  MPa  to  fit  to  the 
experimental  data,  2)  q  =  1.0;  m  =  0.012  -  a  lower  latent  hardening  material,  for  which 
single  slip  behavior  was  found  to  be  given  by  the  Equation  41,  except  for  Ss  which  was 
170  MPa  to  fit  to  the  experimental  data.  The  corresponding  simple  shear  behavior  for  these 
two  cases  are  given  in  Figure  12  with  the  behavior  from  single  slip  define  by  the  Equation 
41  and  experimental  data. 


COMPARISON  OF  EXPERIMENTAL  RESULTS  WITH  MODEL 

PREDICTIONS 


The  suitable  selection  of  the  parameters  of  the  single  slip  behavior,  with  m  =  0.012 
and  q  =1.4,  made  it  possible  to  obtain  a  simulated  uniaxi^  compression  stress-strain  curve 
that  closely  follows  the  experimental  results  as  given  in  Figure  2.  The  parameters  chosen 
to  fit  the  compression  data  remain  fixed  for  the  torsion  simulations.  The  predicted  stress- 
strain  curve  is  dependent  upon  the  rate  sensitivity  as  well  as  the  amount  of  latent  hardening: 
the  higher  the  rate  sensitivity,  the  lower  the  stress  at  a  give  strain  value  and,  the  lower  the 
amount  of  latent  hardening,  the  lower  the  stress  for  a  given  strain  value. 

The  stress-shear  strain  curve  from  the  simulation  of  simple  shear  is  given  in  Figure 
10  with  the  fixed-end  torsion  stress-strain  data.  In  both  simulation  and  experiment,  the 
axial  stresses  are  compressive,  thus  qualitatively  both  have  the  same  results;  but  the 
simulation  over  predicts  the  measured  axial  stress.  Also,  the  experimental  axial  stress 
becomes  tensile  at  higher  shear  strains.  The  simulation  does  not  show  this  tiend  up  to  the 
maximum  simulated  shear  strain  of  2,8.  The  experimental  shear  stresses  are  higher  than 
the  simulated  values  at  low  and  high  strains,  but  the  simulated  stress  is  higher  than  the 
experimental  values  for  intermediate  strain  levels.  This  observation  is  different  from  that 
observed  by  Bronkhorst  et  al.'^  They  observed  a  better  correlation  between  results  from 
experiments  and  simulation.  Also  we  should  note  that  the  two  materials  are  different  since 
the  annealing  temperature  was  different  in  both  studies.  An  annealing  temperature  of  850°C 
was  used  in  their  studies. 

The  other  difference  between  the  two  studies  is  the  annealing  point  in  the  study: 
i.e.,  whether  the  annealing  was  performed  before  or  after  machining  of  the  specimen.  In 
Bronkhorst  study,  specimens  have  been  annealed  after  the  machining,  but  in  this  study  the 
specimens  were  machined  after  annealing  so  as  not  to  distort  the  thin  wall  of  the  specimen. 
With  the  assumption  that  the  specimen  may  have  been  work  hardened  during  machining, 
the  following  two  approaches  were  taken:  1)  a  simulation  was  conducted  by  increasing  the 
initial  single  slip  resistance,  Sq,  so  that  the  initial  portion  of  the  predicted  shear  stress-shear 
strain  curve  matches  the  experimental  results,  and  2)  one  specimen  was  annealed  before 
testing.  Figure  13  gives  a  comparison  of  the  simulated  results  from  the  first  case  with  the 
experiment,  where  the  given  simulation  was  conducted  with  the  initial  single  slip 
resistance,  Sq  =  45  MPa.  The  discrepancy  between  the  experiment  and  the  simulation  at  the 
lower  strains  disappeared  with  this  simulation,  but  the  differences  at  the  intermediate  and 
large  strains  remain  unchanged.  In  the  second  case,  where  the  specimen  was  annealed 
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before  testing,  the  experimental  results  are  given  with  the  simulated  results  in  Figure  14.  In 
this  case  also,  the  differences  in  shear  stress  at  lower  strains  disappeared;  however,  the 
differences  between  the  experiment  and  the  simulation  increased  in  the  intermediate  strain 
range.  At  the  range  of  large  strains,  the  experimental  shear  stress  seems  to  asymptotically 
approach  the  simulated  results.  The  new  simulated  and  experimental  axial  responses  did 
not  appreciably  deviate  from  the  original  simulated  and  experimental  axial  responses, 
respectively. 

Figure  1 1  shows  a  comparison  of  the  simulated  shear  and  axial  stress- shear  strain 
responses  for  an  increase  of  rate  sensitivity  (m  =  0.012  ->  0.04)  and  for  a  decrease  of  latent 
hardening  (q  =  1.4  ->  1.0)  parameters.  The  shear  stress  response  did  not  change 
appreciably,  but  there  is  a  change  in  the  axial  and  hoop  stress  responses.  As  the  strain  rate 
sensitivity  was  increased,  the  axial  and  hoop  stresses  decreas^  for  a  given  strain.  For 
strains  less  than  0.6,  the  axial  and  hoop  stresses  are  higher  with  reduced  latent  hardening, 
but  are  lower  for  reduced  latent  hardening  for  strains  greater  than  0.6. 

In  contrast  to  the  simulations  given  in  Figure  11,  Figure  12  gives  the  simulated 
stress- strain  response  for  both  higher  rate  sensitive  (m  =  0.04)  and  lower  latent  hardening 
(q  =  1 .0)  materials  by  obtaining  the  single  crystal  slip  parameters  for  both  cases  from  the 
compression  data  as  described  in  the  previous  section.  These  simulations  further  deviate 
from  the  experimental  results,  therefore  indicating  that  the  best  simulation  is  obtained  for  m 
=  0.012  (the  experimental  m  from  polycrystal  jump  test  by  Bronkhorst  et  al.)  and  q  =  1.4 
(the  value  suggested  by  Kocks  et.al.'^  and  Asaro  et  al.^  after  reviewing  a  great  number  of 
single  crystal  experiments). 

The  discrepancies  between  predicted  and  experimental  stress-strain  responses  may 
be  due  to  the  followings  reasons: 

1) .  The  Taylor  model  does  not  take  into  account  the  interaction  between  grains.  Also  the 
stress  equilibrium  stress  between  the  grains  is  violated. 

2) .  The  fixed-end  torsion  tests  do  not  represent  simple  shear  perfectly;  this  could  explain 
the  reduction  in  the  magnitude  of  the  experimental  axial  stress.  Lipkin  et  al.'^  have 
experimentally  observed  the  presence  of  axial  strains  in  the  gage  section,  thus  reducing  the 
magnitude  of  the  axial  stress  due  to  partial  relaxing  of  the  axial  constraints.  Also  finite 
element  calculations  by  Lipkin  et  al.'*  and  White^^  show  the  presence  of  axial  strains  in  the 
gage  section.  These  finite  element  calculations  also  show  that  the  actual  shear  strain  in  the 
gage  section  is  lower  than  that  computed  from  the  relative  rotation  measured  at  the  testing 
machine  (grips).  In  addition,  the  finite  element  calculations  show  that  even  for  thin  wall 
specimens,  stresses  are  not  uniform  through  the  thickness  or  close  to  the  shoulder  section 
of  the  gage  area.  Also,  we  should  note  that  these  finite  element  calculations  have  been 
conducted  with  assumed  constitutive  models. 


1 5.  KCXTKS,  U.  F.  The  Relation  Between  Polycrystal  Deformation  and  Single-Ci  ystal  Deformation.  Mel. 
Trans.,  v.  1,  1970,  p.  1121-1143. 

16.  LIPKIN.  J..  CHIESA,  M.  L.  AND  BAMMANN,  D.  J.  Thermal  Softening  of  304L  Stainless  Steel: 
Experimental  Results  and  Numerical  Simulations.  Proc.  of  IMPACT’87,  Bremen,  FRO,  1987. 

17.  WHITE,  C.  S.  An  Analysis  of  the  Thin-Walled  Torsion  Specimen,  submitted  to  the  ASME  J.  Eng. 
Mat  Tech..  1991. 
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STRESS  (MPa) 


Figure  12.  Simulated  Simple  Shear  Stress-Strain  as  a  function  of 
Strain  Rate  Sensitivity  (m  =  0.012  and  0.04)  and 
Latent  Hardening  (q  =  1.4  and  1.0) 

(each  single  slip  behavior  was  deduced  from  compression  experimental  data) 
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STRESS  (MPa) 


Figure  13.  Simulated  Stress-Shear  Strain  with  sq  =  45  MPa  versus 
Experimental  Fixed-end  Torsion  Stress-Shear  Strain  Data 
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STRESS  (MPa) 


SHEAR  STRAIN 


Figure  14.  Simulated  Stress-Shear  Strain  versus  Experimental  Fixed-end 

Torsion  Stress-Shear  Strain  Data 
(specimen  was  heat  treated  after  machining) 
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PREDICTED  TEXTURES 


As  described  in  the  section  entitled  The  Taylor-type  Polycrystal  Model,  the 
evolution  of  the  texture  can  be  obtained  from  the  polycryst^  simulation.  The  initially 
assumed  random  distribution  of  the  crystals  are  given  in  Figure  7  as  equal-area  projection 
{111}  pole  figures.  Predicted  pole  figures  for  both  uniaxial  compression  and  simple  shear 
experiments  are  discussed  and  compared  with  the  corresponding  ideal  textures  in  Figures 
15  and  16  in  the  following  sections.  Experimental  pole  figures  are  being  obtained  and  will 
be  discussed  and  compart  with  the  predicted  pole  figures  in  a  future  report. 

Uniaxial  Compression 

At  a  uniaxial  compression  strain  of  1.3,  equal-area  projection  (111),  (110)  and 
( 100}  pole  figures  of  the  polycrystal  are  given  in  Figure  17.  From  these  pole  figures,  it 
can  be  seen  that  the  polycrystal  has  developed  a  texture  which  is  symmetric  with  the 
loading  axis.  As  seen  by  others  (Bronkhorst  et  al.^'*)  in  uniaxial  compression  tests,  the 
(110}  planes  have  aligned  perpendicular  to  the  loading  axis  (see  the  (110}  pole  in  the 
{110}  pole  figure).  The  predicted  uniaxial  compression  pole  figure  is  identic^  to  the  ideal 
( 1 10}<wvvt'>  fibre  texture  given  in  Figure  15,  where  { 1 10}  planes  align  perpendicular  to 
the  loading  axis  (es-direction). 

Fixed-end  Torsion  (Simple  Shear) 

Figure  18  gives  the  predicted  (111).  1110}  and  {100}  pole  figures  for  a  simple 
shear  test  at  a  shear  strain  of  2.8.  These  pole  figures  are  similar  to  the  ones  predicted  by 
Harren  et  al.*  using  the  time  integration  algorithm  by  Asaro  and  Needleman. 

Ideal  textures  that  are  formed  in  shear  deformation  have  been  classified  in  detail  by 
Canova  et  al.**,  Montheillet  et  al.^^  and  Harren  et  al.®  Following  the  notation  by  Canova  et 
al.,  the  ideal  textures,  [hkl]<uvw>,  where  {hkl}  planes  aligning  with  the  global  shear 
plane  (e2-direction)  and  <uvw>  directions  aligning  with  the  global  shear  direction  (ei- 
direction),  that  form  under  shear  can  be  divided  into  three  categories:  A  fibre  texture  - 
{ 1 1 1  }<«vw>,  B  fibre  texture  -  [hkl]<l  10>,  and  C  texture  -  {001  )<1 10>.  These  ideal 
shear  textures  are  given  as  {111}  pole  figures  in  Figure  18.  It  has  been  hypothesized  by 
Harren  et  al.  that  the  B  fibre  and  the  C  orientation  textures  do  not  produce  axial  stresses, 
but  the  A  fibre  texture  produces  either  compressive  or  tensile  axial  stresses  depending  on 
whether  the  A  fibre  is  either  A'  partial  fibre  or  A'*'  partial  fibre  (partial  A  fibres  are  identified 
in  Figure  16a),  respectively.  Also  in  the  Figure  16(b),  B  fibre  is  divided  into  two  partials. 
By  and  B". 

Comparing  the  predicted  and  ideal  {111}  pole  figures,  prediction  shows  the 
presence  of  a  strong  A'  fibre  and  the  absence  of  A+  fibre.  Also  the  By  fibre  and  the  C 
orientation  are  strongly  present  in  the  pole  figure;  however,  B"  fibre  is  absent  in  the 
prediction. 


18.  CANOVA,  G.  R.,  KOCKS,  U.  F.  AND  JONAS,  J.  J.  Theory  of  Torsion  Texture  Development.  Acta 
Mctall.,  V.  32,  1984,  p.  211-226. 

19.  MONTHEILLET,  F,  COHEN,  M,  AND  JONAS,  J.  J.  Axial  Stresses  and  Texture  Development 
During  the  Torsion  Testing  of  Al,  Cu  and  a-Fe.  Acta  Meiall.,  v.  32,  1984,  p.  2077-2089. 
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Ideal  Textuie 


X2 


Ideal  Textures  in  Uniaxial  Compression 
(110)  <uvw>  Fibre 


Figure  15.  Ideal  Texture  {111}  Pole  Figure  for  Uniaxial  Compression 


X2  X2 


a)  A  Fibre  {lll)<uvw>  -  Ideal  Textures  in  Shear  b)  B  Fibre  (hkl)<II0>  and  C  {001)<110>  -  Ideal  Textures  in  Shear 
(A-  :  solid  lines.  A-*-  :  dashed  lines  and  A  :  dots)  (By  :  solid  lines.  Bn  :  dashed  lines  and  C  :  dots) 


a).  A  Fibre 

(A-*^  :  Solid  Lines  and  A*  :  Dashed 
Lines) 


b).  B  Fibre  and  C 
(By  :  Solid  Lines,  B"  :  Dashed 
Lines  and  C  :  Dots) 


Figure  16.  Ideal  Texture  (111)  Pole  Figures  for  Simple  Shear 
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□ 


(111)  Pol«  flgarc  la  aaiailal  caafrcMloa  at  ttraia 


X2 


1.3 


D  (110)  Pol<  flgare  (a  aaiailal  conprctiioa  at  itrala  s  1.3 

(■a0.0n;4scl.4) 

X2 


□  (100)  Pale  rigvre  la  aaiailal  coaiprcfaloa  at  straia  s  1.3 

(m=0.0I2;<|=I.4) 


Figure  17.  Computed  Pole  Figures  from  Rate  Dependent  Taylor  Model  for 

Uniaxial  Compression 
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The  presence  of  A‘  fibre  is  responsible  for  the  compressive  axial  stresses  seen  in  the 
predicted  as  well  experimental  fixed-end  torsion  tests.  The  perturbations  in  the 
compressive  axial  stress  with  increasing  strain  is  due  to  the  migration  of  orientations 
fronVto  A-  to/from  By  fibre  and  C  orientation.* 

SUMMARY 

Large  deformation  uniaxial  compression  and  fixed-end  and  ffee-end  torsion  (simple 
shear)  experiments  were  conducted  for  annealed  OFHC  Copper  to  obtain  stress-strain 
behavior  of  the  material.  The  stress-strain  behavior  for  these  experiments  were  also 
predicted  using  a  Taylor  type  rate  dependent  polycrystal  model  by  Asaro  and  Needleman^. 
Simulation  of  these  experiments  was  conducted  using  a  recently  developed,  highly 
efficient,  fully  implicit  time  integration  scheme  by  Kalidindi  et  al.^  In  the  initial  phase  of 
the  simulation,  the  evolution  of  the  constituent  single  crystal  slip  system  deformation 
resistance  was  estimated  using  the  experimentally  determined  compressive  stress-strain 
behavior  of  the  polycrystal.  With  this  crystal  constitutive  behavior,  the  stress-strain 
behavior  of  the  polycrystal  for  simple  shear  was  computed.  In  addition,  evolution  of  the 
crystallographic  texture  was  computed  for  both  compression  and  simple  shear  tests. 

Experimental  shear  stress  was  insensitive  to  changing  from  fixed-end  to  ffee-end 
condition  of  testing.  The  predicted  shear  stress-strain  behavior  for  simple  shear  test 
compared  reasonably  well  with  the  experimental  results.  The  simulated  axial  stress 
response  agrees  qualitatively  (compressive)  with  the  experimental  observations,  but  over 
estimates  the  experimentally  observed  axial  compressive  stress. 

In  uniaxial  compression  tests,  the  simulated  stress- strain  behavior  decreases  with  an 
increase  of  strain  rate  sensitivity  or  decrease  of  the  latent  hardening  for  the  same  strain.  In 
simple  shear,  shear  stress  response  is  not  sensitive  to  changes  in  either  strain  rate 
sensitivity  or  latent  hardening.  The  axial  and  hoop  stresses  decrease  with  the  increase  in 
strain  rate  sensitivity.  With  the  decrease  of  latent  hardening,  the  axial  and  hoop  stresses 
increase  at  lower  shear  strains,  but  decrease  at  higher  shear 

The  predicted  texture  of  grains  for  compression  tests  assumes  an  orientation  given 
by  the  ideal  { 1  10}<mvw>  fibre  where  {110}  planes  become  perpendicular  to  the  loading 
axis.  For  simple  shear,  predictions  show  the  presence  of  the  ideal  A*  (( 1 1 1  }<wvw>  panial 
fibre)  and  By  {[hkl}<l  10>  partial  fibre)  fibres  and  C  ((001  )<1 10>)  orientation,  where  in 
this  notation,  planes  are  parallel  to  the  shear  plane  and  directions  are  along  the  shear  - 
direction. 
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APPENDIX 


Computer  Program  PolyCrystal-MTL  (PC-MTL-Ver:  01). 
A  C  version  of  the  program  will  be  available  in  future. 

program  Tpoiycrytal; 

{••  written  by  Tusit  Weerasooriya  of  Army  Material  Technology  Lab  **) 

{*•  December  1990  **} 

(This  program  computes  the  stress  history  for  a  given  history  (including  ) 

{reverse  loading)  of  deformation  gradient  tensor  (F).  Continuum  is  ) 

(considered  as  a  poly-crystal  of  n  number  of  crystals  with  initially  known  ] 
(orientation  (ex.  randomly).  The  program  also  tracks  the  evolution  of  ) 
(anisotrophy  (crystalogiaphic-texture)  with  the  deformation.  State  variables) 

{ -  crystal  slip  resistance  and  Fp  -  for  each  crystal  are  tracked  and  are  } 

(saved  in  a  restart  file.  Using  the  restart  Tde,  deformation  can  be  further) 

(extOKlcd) 

const 

ncrystals  =  1(X); 
startcrystal  =  1; 
ntimesteps  =  4(X); 

lastcrystal  =  startcrystal  +  ncrystals  -  1; 
lotcrystals  =  1(X); 

qq=  1.0; 

maxstrain  =  2.80; 
minsuain  =  0.000; 
startstrain  =  0.000; 
sdgama  =  O.OI; 

ho  =  700.0; 
ss  =  155.0; 
so  =  16.0; 
aa  =  3.80; 

gamadoto  =  0.001; 
m  =  0.012; 
straindot  =  0.001; 
mu  =  44(XX).0; 

type 

table66  =  array[1..6,  1..6)  of  extended; 
tablc33  =  array [1.. 3, 1..31  of  extended; 
rtable33  =  array [1.. 3,  1..3)  of  real; 
vector3  =  array  (1.. 3]  of  extended; 
vectoi6  =  array[1..6)  of  extended; 
vectorl2  =  array[1..12)  of  extended; 

(ablel2l2  =  array[1..12,  1..I2]  of  real; 

Scl2  =  array[1..12)  of  table33; 

Qvector  =  array)  1  ..lotcrystals]  of  table33; 
svcctor  =  array) l..totcry.stalsl  of  vcctorI2; 


(mult  for  latent.h  for  noncoplanar) 

(max  Slop  strain) 

(min  stop  strain  for  reverse) 

(starting  strain  before  reverse) 
(maximum  delta _gamafor  iteration] 

{ s_dot=ho(  1  -s/ssY^ ) 


(reference  strain  rate) 
(suain  rate  sensitivity] 
(strain  rale  of  the  lest) 


34 


Tvector  =  array[l..loicrystals]  of  vector^; 
rtstepv  =  array [l..numestcpsl  of  real; 
tstepv  =  array[l..ntimesteps]  of  extended; 
pi  1 1  =  array (1.. 8]  of  vector3: 
pi  10  =  array[1..12]  of  vectorS; 
plOO  =  array [1.. 6]  of  vector3; 

var 

answer  string[3]; 

restart:  string[2]; 

lesttype,  testdim,  strainlimit:  char, 

crystal,  il,  Itimestep,  n,  i,  j,  itemp:  integer; 

qab:  tablel212; 

k,  r3kin2niud6,  r3kd2,  sod3,  temp,  scaler:  real; 

Rt,  cdgamamax,  newdgamamax,  ptsdgamamax,  dgamamax:  real; 
plllc,plllg:  pill; 
pi  10c,  pllOg:  pi  10; 
plOOc,  plOOg:  plOO; 

Sc,  Sa:  Sc  12; 

Q:  table33; 
eulera:  vector3; 

166,  null66:  table66; 

133,  null33,  temp33,  temp  133,  temp233,  Feinter:  table33; 
nullv6:  vector6; 
nullvl2,  sovl2:  vector  12; 

Qptr  'HJveclor, 

Fpiptr,  nFpiptr  XJvector; 
sptr,  nsptr;  ^svector; 

Tsvptr,  nTsvptr  ''Tvector; 

X,  y,  msdot:  real; 
time,  strain:  rtstepv; 

T1  lavptr,  T22avptr,  T33avptr,  T12avptr:  ''tstepv; 
deltat:  real; 

Tfile,  pFile,  pOfile,  pOOfile,  pifile,  piOFile,  piOOfile:  text; 
var 

l,  iterate:  integer; 
alpha,  beta:  integer, 

F,  Fe,  Fpi,  Fpinew,  R,  InU:  lable33; 

A,  Ts,  Tstr,  Tcauchy:  table33; 

Tsnewvec,  Tvecerror,  Tstrvec,  Tfn,  tcmp6:  vcctor6; 
pole,  newpole:  veclor3: 

Ba,  Ca:  Sc  12; 

Savec,  Cavec:  array[1..121  of  veclor6; 

Jacobian,  Uac,  temp66:  tablc66; 

s,  snew,  serror,  st,  iau,  dellagama,  tcmpl2;  vcctorl2; 

gamopdelt,  si,  s2,  s3,  s4,  s5:  extended; 

Tout  sout,  iterouc  boolean; 


procedure  mprint  (am:  tablc33; 

n:  integer); 
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(print  3x3  in  MatrixForm] 
var 

i,  j:  integer, 
begin 

for  i  :=  1  to  n  do  begin 
writeln; 

for  j  :=  1  to  n  do 
write(^[i,  j] :  10  :  3); 
end;  (‘for  i*) 
writeln; 
end; 

{ - ) 

procedure  vprint  (av;  vectorl2; 

n:  integer); 

(print  a  12  element  vector) 
var 

i,  j:  integer; 
begin 

for  i  :=  1  to  n  do 
write(av[i] :  12  :  6); 
writeln; 
end;  (vprint) 

{ - ) 

procedure  m6print  (am:  table66; 

n:  integer); 

(print  a  6x6  in  MatrixForm) 

var 

i,  j:  integer, 
begin 

for  i  :=  1  to  n  do  begin 
writeln; 

for  j  :=  1  to  n  do 
write(am[i,  j] ;  10  :  4); 
end;  (*for  i*) 
writeln; 
end; 

{ - ) 

procedure  v6print  (av:  vector6; 

n:  integer); 

(print  vector  of  6  elements) 
var 

i,  j:  integer; 

begin 

for  i  :=  1  to  n  do  begin 
write(av[i) :  12:6); 
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end;  (*for  i*) 
writeln; 
end;(vprint) 

{ - ) 

procedure  inverse  (am:  table33; 

var  bm:  table33; 
n:  integer); 

(inverse  a  3x3  matrix  by  back  substitution) 
var 

i.  j,  k:  integer, 
p:  extended; 
s;  extended; 

begin 

for  i  :=  1  to  n  do 
for  j  :=  1  to  n  do 
bm[i,  j]  :=am[i,  j); 


for  i  :=  1  to  n  do  begin 
p  :=  bm[i,  i); 

p  :=  1.0  /  p; 
bm[i,  i]  :=  1.0; 
for  j  1  to  n  do 
bm[i,  j]  ;=:bm[i,  j]  •  p; 
for  k  ;=  1  to  n  do 
if  (k  <>  i)  then  begin 
s  ;=  bm[k,  i); 
bm[k,  i]  :®  0.0; 
for  j  :s  1  to  n  do 
bm[k,  j]  :=  bmlk,  j]  -  s  *  bm[i,  j]; 
end;(*ir*) 
end;  (fori) 
end;  ( inverse) 


procedure  inverseb  (am:  tablc66; 

var  bm:  table66; 
n:  integer); 

( inverse  a  6x6  matrix  by  back  substitution ) 
var 

i,  j,  k:  integer; 
p:  extended' 
s:  extended; 

begin 

for  i  :=  1  to  n  do 
for  j  :=  1  to  n  do 
bm[i,  j)  :=  am[i,  j); 
for  i  :=  1  to  n  do  begin 
p  :=  bm(i,  i); 
p  :=  1.0 /p; 
bm[i,  i)  :=  1.0; 
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for  j  :=  1  to  n  do 
bm[i,  j]  :=  bm[i,  j]  ♦  p; 
for  k  :=  1  to  n  do 
if  (k  <>  i)  then  begin 
s  :=  bin[k,  i]; 
bin[k,  i]  :=  0.0; 
for  j  :=  1  to  n  do 
bm[k,  j]  :=  bm[k,  j]  -  s  *  bm[i,  j]; 
end;(*if*) 
end;  (for  i} 
end;{inverse6) 

{ - 

procedure  aikbkj  (am:  table33; 

bm:  table33; 
var  cm:  table33); 

(Aik  Bkj  tensor  muldplicadon] 
var 

i,  j,  k:  integer, 
begin 

for  i  :=  1  to  3  do 
for  j  ;=  1  to  3  do  begin 
cm[i,  j]  :=  0.0; 
for  k  :=  1  to  3  do 

cm[i,  j]  :=  cm(i,  j]  +  am(i,  k]  *  bm[k,  j]; 
end;  (fork] 
end;(aikbkj} 

( - ) 

procedure  aikbjk  (am:  table33; 

bm:  table33; 
var  cm:  tabie33); 

(Aik  Bjk  tensor  muldplicadon  •  A  (B)T} 
var 

i,  j,  k:  integer; 
begin 

for  i  :=  1  to  3  do 
for  j  ;=  1  to  3  do  begin 
cm[i,  j]  :=  0.0; 
for  k  :=  1  to  3  do 

cm[i,  j]  :=  cm[i,  j]  +  am[i,  k]  •  bm[j,  k]; 
end;  (fork) 
end;  (aikbjk) 

{ - ) 

procedure  mplusm  (am;  tablc33; 

bm:  table33; 
var  cm:  table33); 

(AiJ+Eij) 

var 

i,  j:  integer; 
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begin 

for  i  :=  1  to  3  do 
for  j  :=  1  to  3  do 
cm[i,  j]  :=  ain[i,  j]  +  bin[i,  j]; 
end;(mplusm) 

{ - - - ) 

procedure  mdotv  (am:  table33; 

av:  vector3; 
var  cv:  vector3); 

(AijBj) 

var 

i,  k:  integer, 
begin 

for  i  ;=  1  to  3  do  begin 
cv[i]  :=  0.0; 
for  k  :=  1  to  3  do 
cv[i]  :=  cv[i]  +  am[i,  k]  *  av(lc]: 

end;  {fori} 
end;(mdotv} 

( - ) 

procedure  smulm  (scalerl:  real: 

am:  table33; 
var  bm:  table33); 

{s*Aij} 

var 

i,  j:  integer, 
begin 

for  i  :=  1  to  3  do 
for  j  :=  1  to  3  do 
bm[i,  j]  :=  scalerl  ♦  am[i,  j]; 
end;  (smulm) 

{ - ) 

function  absmax  (av:  vectorl2):  real; 

(abs(max(Aij))) 

var 

i:  integer, 
begin 

temp;=abs(av[l]); 

for  i  :=  2  to  12  do 
if  (abs(av[i]))  >  temp  then 
temp;=  abs(av[ij); 

absmax  :=  temp; 
end;  (absmax) 

function  deG3  (am:  table33):  real; 
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var 

i:  integer, 
begin 

deG3  :=  am[l,  1]  *  (am[2. 2]  *  am[3. 3]  -  am[3, 2]  •  am[2. 3])  -  am[2. 1]  *  (am[l.  2]  *  am[3, 3]  -  am[3, 2]  *  am[l. 
3])  +  ain[3. 1]  *  (amil,  2]  •  am[2. 3]  -  am[2. 2]  •  ain[l.  3]); 

end:{(let33) 

{ - j 

procedure  vouterprod  (var  av,  bv:  vector3; 

var  cm:  table33): 

( Ai  Bj  tensor  product) 
var 

i,  j:  integer, 
begin 

for  i  :=  1  to  3  do 
for  j  :=  1  to  3  do 
cm[i,j]  :=av[i]  *  bvlj]; 
end;  {vouterprod) 

{ - ) 

procedure  genSc  (var  Sc:  Sc  12; 

var  pi  11c;  pill); 


var 

sqrt2,  sqrt3:  real; 
i,  j:  integer, 

sll,  sl2,  sl3,  s21,  s22,  s23,  s31,  s32,  s33,  s41,  s42,  s43:  vector3; 
begin 

sqrt3  :=  sqrt(3.0); 
sqrt2  ;=  sqit(2.0); 
for  i  ;=  1  to  3  do  begin 
plllc[l,  i]  :=  1.0/sqrt3;  (1,1,1) 

plllc[2,  i]  :=  1.0/sqrt3; 
plllc[3,  i]  :=  1.0/sqrt3; 
plllc[4,  i]  :=  1.0/  s<irt3; 


end; 

plllc[2,l]:=  -plllc(2, 1); 

(-1.1,1) 

plllc(3, 2]  :=  -plllc[3, 2); 

{  l.-l.l) 

plllc[4, 3]  :=  -plllc[4, 3); 
for  i  :=  1  to  3  do  begin 

{  l.l.-l) 

plllc(5,  i]  :=  -plllc[l,i]; 

(-1.-1.-1) 

plllc[6,  i]  :=  -plllc[2,  i); 

{  l.-l.-l) 

plllc[7,  i]  :=  -plllc[3,  i); 

(-l.l.-l) 

plllc(8,  i]:=  -pmc(4,i]; 

end; 

slid)  := -1.0 /sqrt2; 
sll[2J  :=  1.0/sqrt2; 

(-1.-1. 1) 

si  1(3]  :=  0.0; 
sl2(l]  :=-1.0/sqrt2; 
sl2(2)  :=  0.0; 

(-1.1.0) 
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sl2[3]  :=  1.0/sqrt2; 

1-1,0. 1) 

sl3[l]  :=  0.0; 
sl3[2]  :=-1.0/sqrt2; 
sl3[3]  :=  1.0/sqrt2; 

(O.-l.l) 

s21[l]  :=1.0/s(pt2; 
s21[2]  :=1.0/sqrt2; 
s21[3]  :=  0.0; 

{  1.1.0} 

s22[l];=1.0/sqrt2; 
s22[2]  :=  0.0; 
s22[3]  :=  1.0/sqrt2; 

{  1.0.1) 

s23ll]  :=  0.0; 
s23[2]  ;=.1.0/sqrt2; 
s23[3]  :=  1.0/sqrt2; 

(O.-l.l] 

s31[l]  :=  1.0/sqrt2; 
s31[2]  ;=  l.O  /  sqrt2; 
s31[3]  :=  0.0; 

{  1,1,0) 

s32[l]  :=  0.0; 
s32[2]  ;=  1.0/sqn2; 
s32[3]  :=  1.0/sqrt2; 

{0.1,1} 

s33[l]  :=  1.0/sqrt2; 
s33[2]  :=  0.0; 
s33[3]  :=  -1.0/sqrt2; 

{  1.0,-1) 

s41[l]  :=  0.0; 
s41[2]  :=  1.0/sqrt2; 

S4U31  ;=  l.O  /  sqrt2; 

{  0.1.1] 

s42[l]  :=  1.0/sqrt2; 
s42[2]  :=  0.0; 
s42l3]  :=  1.0/sqrt2; 

(  1.0.1) 

s43(l]  :=  1.0/sqrt2; 
s43[2]  :=-1.0/sqil2; 
s43[3]  :=  0.0; 

{  l.-l.O) 

for  i  :=  1  to  3  do 
for  j  ;=  1  to  3  do  begin 
Sc[l.i.j]  :=sll[i]*plllc[l,j]: 
Sc[2,  i.j];=sl2m*plHc[l,j]; 
Scp.i.j]  :=sl3[i]*plllc[l,j]; 
Sc[4,  i.j]  :=s21[i]*plllc[2,j]; 
Sc[5,  i,j]:=  s22[i]*plllc[2,j]: 
Sc[6.  i,j]  :=  s23[i]*plllc[2,j]; 
Sc[7.i,j]  :=  s3Un*plllc[3.j]; 
Sc[8.  ij]  :=  s32[i]*plllc[3J]; 
Sc[9, i.j]  :=  s33[i]*plllc[3.j]: 
Sc(10.  i.j]  ;=s41(i]  *pnic[4.j]; 
Sc[ll.i.j]  :=  s42[i]*plllc[4.j]; 
Sc[12.  i.j]  :=  s43[i]*plllc[4,j]; 
end;  ( for  j] 
end;(genSc) 


procedure  genQ  (var  Q:  table33; 

eulerangle:  vector3); 


var 

phi,  omega,  ihrL,:  sau, 

cphi,  sphi,  comega,  somega,  ctheta,  stheta:  real; 
begin 

phi  :=eulerangle[l};  {0=<phi=<360) 

theta  :=  eukrangle[2];  {0=<theta=<180} 

omega  :=  eulerangle[3];  {(>=<omega=<360} 

cphi  ;=  cos(phi); 

sphi  ;=  sin(phi); 

ctheta  ;=  cos(theta); 

stheta  :=  sin(theta); 

comega  :=  cos(omega); 

somega  :=  sin(omega); 

Q[l,  1]  :=  cphi  *  comega  -  sphi  *  somega  *  ctheta; 

Q[l,  2]  :=  sphi  •  comega  +  cphi  *  somega  *  ctheta; 

Q(l,  3]  :=  somega  *  stheta; 

Q(2, 1]  :=  -cphi  *  somega  -  sphi  ♦  comega  *  ctheta; 

Q[2, 2]  :=  -sphi  ♦  somega  +  cphi  ♦  comega  *  ctheta; 

Q[2, 3]  :=  comega  *  stheta; 

Q[3, 1]  :=  sphi  *  stheta; 

Q[3,  2]  :=  -cphi  •  stheta; 

Q[3, 3]  ;=  ctheta; 
end;  {genQ) 

{ - ) 

procedure  eulenead  (totcrystals:  integer); 
var 

i:  integer, 
eafile:  text; 

begin 

open(earile,  'text^OO"); 
reset(eafile); 

temp  :=  3.14159265  /  180.0; 
for  i  :=  1  to  totcrystals  do  begin 
readln(eanie,  euleta[2],  eulera{l],  eulera[3]); 
eulera)!]  :=  temp  *  eulera[l];  (phi) 
eulera[2]  :=  temp  eulera[2];  (theta) 
euleia[3)  :=  temp  *  eulera{3);  (om^) 

genQ(Qptr^[i),  eulcra); 

end;  (fori) 
closefeafile); 
end;  (eulenead) 

{ - ) 

function  trace  (am:  table33):  real; 
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(calculate  the  trace  of  a  matrix  -  Akk) 
var 

i:  integer, 
begin 
scaler  :=  0.0; 
for  i  :=  1  to  3  do 
scaler  :=  scaler  am[i,  i]; 
trace  :=  scaler, 
end;  (trace) 

{ - ) 

procedure  symmtov  (am:  table33; 

var  avec:  vectoi6); 

(vectorise  a  symmetric  matrix) 
var 

i:  integer, 
begin 

for  i  :=  1  to  3  do 
avec(i]  am[i,  i); 
avec[4]  :=am[l,2]: 
avec[5]  ;=am[l,3]; 
avec[6]  :=  am[2, 3); 
end;  (symmtov) 

{ - ) 

procedure  symmtopv  (am:  tabie33; 

var  avec:  vector6); 

(vectorise  a  matrix) 
var 

i,  j:  integen 
begin 

for  i  :=  1  to  3  do 
for  j  :=  1  to  3  do 

temp33[i,  j]  :=  0.5  •  (am(i,  j]  +  am[j,  i]); 


for  i  :=  1  to  3  do 
avec[i]  :=  temp33[i,  i); 
avec(4]  :=  2.0  *  temp33[l,  2] 
avec[5]  :=  2.0  •  temp33[l,  3] 
avec(61  :=  2.0  •  temp33(2, 3] 
end;  (symmtopv) 


procedure  vtosymm  (avec:  vectorb; 

var  »  ‘'le33); 

( convert  the  6  elemeni  ecte  o  the  corresponding  symmetric  3x3  matrix ) 

var 

i,  j:  integer, 

begin 


for  i  :=  1  to  3  do 
am(i,  i]  :=  avec(i]; 
am(l,2]  :=  avec[4]; 


am[l,3]  :=avec[5]; 
am[2, 3]  :=  avec[6]: 
ain[2, 1]  :=  ain[l,  2]; 
ain[3, 1]  :=  ain[l,  3]; 
ainr3  2]  :=  ain[2,  31- 
end;  (vtosymm) 


procedure  eapolexy  (hkl:  vector3;  ' 
var  X,  y:  real); 

(*  equal  area  projection  pole  figure  generation  *) 
begin 

(equal  area-projection  pole  Rgure] 
temp  :=  sqrt(1.0  +  hkl[3]); 

X  :=  hkl[l]  /  temp; 
y  :=  hkl[2]  /  temp; 

(stereographic-projection  pole  figure] 
{temp:=l.(>+hkl[3];} 

{x:=hkl[l]/(l.(>fhkl[3]);) 

{y:=hkl[2]/(l.(>4-hkl[3]);} 

end;  (eapolexy) 

{ - 

procedure  RUdecomp  (F:  table33; 

var  R,  InU:  table33); 

(polar  decomposition  of  F  s  RU) 
var 

c,  c2:  table33; 
i,  j,  k:  integer, 

ic,  iic,  iiic,  iu,  iiu,  iiiu,  ic2,  ic3,  iic2,  iunew:  extended; 
aO,  al.  a2,  a3,  a4,  iu3,  srtiiic:  extended; 
begin 
c  :=  null33; 
for  i  :=  1  to  3  do 
for  j  :=  i  to  3  do 
for  k  :=  1  to  3  do 
c(i,j]  :=c[i,j]  +  F[k,i]*F[k.j]; 
c[2.  l]:=c[l,2]; 
c[3.  l]:=c[l,3]; 
c[3. 2]  :=c(2, 3]; 

(c^l 

c2  :=  null33; 
for  i  ;=  1  to  3  do 
forj:=ito3do 
for  k  :=  1  to  3  do 
c2[i,  j]  :=  c2(i,  j]  +  c[i,  k]  *  c[k,  jj; 
c2[2,  l]:=c2[1.2]; 
c2[3. 1]:=c2[l,3]; 
c2[3.  2]  :=  c2[2, 3]; 

(calculate  inavriants  of  c) 
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ic  :=  c[l,  1]  +c[2, 2]  +  c[3, 3]; 

iic  :=  (ic  •  ic  -  (c2[l,  1]  +  c2[2, 2]  +  c2[3, 3]))  /  2.0; 

iiic  ;=  det33(c); 

{calculate  iiiu=sqit(iiic)) 

(calculate  iu) 
iunew  :=  ic; 
repeat 
iu  :=  iunew; 
iu3  :=  iu  *  iu  *  iu; 
sitiiic  :=  sqrt(iiic); 

iunew  :=  iu  -  (iu  *  iu3  -  2.0  *  ic  *  iu  *  iu  -  8.0  *  iu  *  srtiiic  f  ic  *  ic  -  4.0  *  iic)  /  (4.0  *  iu3  -  4.0  *  ic  *  iu  -  8.0  * 
srtiiic); 

until  abs(iunew  -  iu)  <  0.(X)0001; 
iu  :=  iunew; 

iiu  :=  sqrt(iic  +  2.0  *  iu  *  sqrt(iiic)); 

(lAl) 

for  i  ;=  1  to  3  do 
for  j  :=  1  to  3  do  begin 
temp33[i,  j]  ;=  c[i,  j]  +  (iiiu  /  iu)  *  I33[i,  j]; 
temp233[i,  j]  :=  c[i,  j]  +  iiu  •  I33[i,  j]; 
end;  { for  j} 

inverse(teinp33.  tenipl33. 3); 
aikbkj(templ33,  tenip233,  InU); 
for  i  1  to  3  do 
for  j  ;s  1  to  3  do 
InU(i,j]  :=  InU[i.j]/iu; 
aikbkj(F.  InU.  R); 

{  mprint(R,3);} 

end;(RUdecomp) 

{ - ) 

procedure  initialize; 
var 

crystal,  i,  j:  integer, 
sqrt2:  real; 

rsFpiptn  file  of  Qvector; 
rsTsvptr:  file  of  Tvector; 
rssptn  file  of  svector; 
rsreal;  file  of  real; 

begin 

for  i  :=  1  to  3  do 
for  j  :=  1  to  3  do 
if  (i  =  j)  then  begin 
I33[i.j]  :=  1.0: 
nuU33(i,  j]  :=  0.0; 
end 

else  begin 

I33[i.j]  ;=  0.0: 
null33[i,  j]  :=  0.0; 


45 


end; 

for  i  :=  1  to  6  do 
for  j  :=  I  to  6  do 
if  (i  -  j)  then  begin 
I66fi,  j]  :=  1.0: 
niill66[i,  j]  :=  0.0; 
end 

else  begin 

I66Ii,  j]  :=  0.0; 
null66[i,  j]  :=  0.0; 
end;  (else) 
for  i  :=  1  to  6  do 
nullv6[i]  :=  0.0; 
for  i  :=  1  to  12  do  begin 
nullvl2[i]  :=  0.0; 
sovl2[i]  :=  so; 
end:{i) 

sqrt2  :=  sqtt(2.0); 
for  j  :=  1  to  3  do 
for  i  ;=  1  to  3  do 
pll0c[j,  ij 1  0  /  sqn2; 


pll0c[l.  1]  :=  0.0; 

{0,1,1) 

pll0c(2,2]  0.0; 

(  1,0,1) 

pllOcf3.3]  :=  0.0; 

{  1.1.0) 

for  j  :=4  to  6  do 

for  i  :=  1  to  3  do 

pllOclj,  i]  :=pll0c|j  -  3,  i]; 
pll0c[4, 2]  :=  -pll0c[4, 2]; 

(0,-l,l) 

pll0c[5,  l]:=-pll0c[5. 1]; 

(-1.0.1) 

pll0c[6, 1];=  -pll0c[6. 1]; 

(-1.1.0) 

for  j  :=  7  to  12  do 

(O.-l.-l) 
(-1.0-1) 
(-1.-1. 0) 

for  i  :=  1  to  3  do 

(O.l.-l) 

{  l.o.-l) 

{  l.-l.O) 

pi  lOcfj,  i]  :=  -pi  lOcU  -  6,  i]; 


for  j  ;=  1  to  6  do 
for  i  :=  1  to  3  do 


pl00c[j,  i]  ;=  0.0; 


plOOcd,  1]  :=  1.0; 

{  1.0.0} 

pl00c[2. 2]  :=  1.0; 

{  0.1.0} 

pl00c(3,  3]  :=  1.0; 

(0.0. 1) 

pl00c(4,  1)  :=-1.0; 

(-1.0,0) 

pl00c[5.2]  :*-1.0: 

(O.-l.O) 

pl00c(6, 3]  :=-1.0; 

(O.O.-l) 

for  i  :=  1  to  12  do 
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for  j  :=  1  to  12  do 

if  (((i  >=  1)  and  (i  <=  3))  and  (0  >=  1)  and  (j  <=  3)))  or  (((i  >=  4)  and  (i  <=  6))  and  ((j  >=  4)  and  (j  <=  6)))  or 
(((i  >=  7)  and  (i  <=  9))  and  ((j  >=  7)  and  (j  <=  9)))  or  (((i  >=  10)  and  (i  <=  12))  and  ((j  >=  10)  and  (j  <=  12))) 

then 

qab[i,  j] 1.0  *  he 
else 

tiab[i,j]  :=qq*ho; 


if  (restart  =  'nn*)  then  begin 

open(piflle,  'pil  1 1 

rewhte(pifile); 

open(piOrile,  'pill0.pr2'); 

rewiite(piOfile); 

open(pi00nie,  'pil00.pr2'): 

rewrite(pi00file); 

rewriteCTfile): 

rewrite(pfile); 

rewrite(pOfile); 

rewrite(p00filc); 

time[l]  :=  0.0; 
strain[l]  :=  0.0; 

T12avptr^[l];=0.0; 

Tllavpti^fl]  :=  0.0; 

T22avptr^(l]  :=  0.0; 

T33avptr^(l]  :s0.0; 

writelnCmie, '  STRAIN  T12av  T33av  T22av  Tllav  T33av-p ); 

writeln(Tfile,strain[l] :  12 : 4.T12avptr^[l] ;  12  :  4,T33avptr^[l] :  12 :  4,  T22avptr''[l] :  12  :  4,  T1  lavptr''[l] :  12 
:  4,  -(T33avptr^[l]  -  (T1  lavptr''[l]  +  T22avptr''[l])  /  2.0) :  12  :  4); 

for  crystal  :=  startcrystal  to  lastcrystal  do  begin 
sptr''[crystal]  :=  sovl2; 

Fpiptr^[crystal]  :=  133; 

Tsvptr''(crystal]  :=  nullv6; 

writeln(pirile, "); 
for  i  ;=  1  to  8  do  begin 
indotv(Qptr''[crysial],  pi  1  lc[i],  pi  1  lg[i]); 
if  plllg[i,  3]  >=  0.0  then  begin 
eapolexy(plllg(i],x,  y); 
write(pifile,  x  :  10  :  4,  y  :  10  :  4); 
end;  (if) 
end;  (fori) 

wrileln(pi0rile,  ”); 
for  i  :=  I  to  12  do  begin 
mdolv(Qptr''(crystal],  pi  10c[i],  pi  lOgli]); 
if  pllOgfi,  3]  >=  0.0  then  begin 
eapolexy(pll0g(i],x,  y); 
write(pi0file,  x  :  10  :  4,  y  ;  10  :  4); 
end;  (if) 
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end;  (fori) 

writeln(piOOfile, "); 
for  i  ;=  1  to  6  do  begin 
mdotv(Qptr*  [crystal],  pl00c(i],  pl00g[il); 
if  pl00g[i,  3]  >=  0.0  then  begin 
eapolcxy(pl00g[i].  x,  y); 
write(pi00file,  x  :  10  :  4,  y  :  10  :  4); 
aid;{if) 
end;  (fori) 
end;  (crystal] 
close(pinie); 
close(piOrile); 
close(pi00file); 

end  (if) 

else  if  (restart  =  ’yy")  or  (restan  =  'ry')  then  begin 

open(rsFpiptr,  'rsFpptr'): 
open(rssptr,  'rssptr'); 
opcn(rsTsvptr,  'rsTsvptr'); 
open(rsreal.  'rsreal'); 

reset(rsFpiptr); 

reset(rssptr); 

reset(rsTsvptr); 

reset(rsreal); 

read(rsFpiptr,  Fpiptr''); 
read(rssptr,  sptr'); 
read(rsTsvptr,  Tsvptr^); 
read(rsreal,  time[l],  strain[I],  deltat); 

close(rsFpiptr); 

close(rssptr); 

close(rsTsvptr); 

close(rsreal); 

{  append(Tfile);  for  PC) 

seekCrrUe,  maxlongint);  ( for  Mac  ] 

rewrite(pfile); 
rewrite(pOfile); 
rewrite(p(X)file); 
end;  (else) 

end;  (initialize) 

{ - 

procedure  saversfiles; 

(save  intermediate  data  for  restarting  the  calculation) 
var 

rsFpiptr:  file  of  Qvector; 
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rsTsvptr:  file  of  Tvector; 
rssptr  file  of  svector; 
rsreal:  file  of  real; 

begin 

open(rsFpiplr,  'rsFpptr'): 
open(rssptr,  'rssptr^; 
open(rsTsvptr,  'rsTsvptr'); 
open(rsreal,  'rsreal'); 

rewrite(rsFpiptr); 

iewTite(issptr); 

rewTite(rsTsvptr); 

iewrite<rsieal); 

write(rsFpiptr,  Fpiptr^'); 
write(rssptr,  sptr^); 
write(rsTsvptr,  Tsvptr^); 

wnte(rsreal,  time[ltimestep],  strain[ltimestcpl,  deltat); 

close(rsFpiptr); 

close(rs^tr); 

close(rsTsvptr); 

cIose<rsieal); 

end;{saversriles} 

{ - ) 

function  sign  (etemp:  extended):  real; 

(if  negative  return  -1.0  and  otherwise  return  I.O) 
var 

si:  real; 
begin 

if  etemp  <  0.0  then 
si  :*-1.0 
else 

si  :=  1.0; 
sign  :=:  si; 
end;  (sign) 


procedure  tau.dgama; 

(calculate  tau  and  delta_gama  for  each  slip  system } 
var 

i,  alpha:  integer, 

si,  s2,  s3,  s4,  s5:  extended; 

begin 

tau  :=  nullvl2; 

for  alpha  :=  1  to  12  do  begin 
for  i  :=  I  to  6  do 

tau(alpha]  :=  tau(alpha]  +  nTsvptr^[crystal](i]  *  Savcc{alpha,  i]; 
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if  (abs(tau[alpha])  <  so  /  100000.0)  then 
deltagamafalpha]  :=  0.0 

else  begin 

s3  :=  abs(Lui[alphal  /  nspti^[crystal,  alpha]); 
s4  :=  In(^)  /  m; 

deltagamaCalpha]  gamopdelt  *  sign(tau[alpha])  *  exp(s4); 
end;  (if) 
end;  (for  alpha) 

end;(tau_dgania} 

{ - } 


label 

10,  20; 

{**'''**********endof  procedures,  functions  and  labels************) 


begin 

k  :=  3.11S  *  mu;  (k- bulk  modulus) 

r3km2mud6  :=  (3.0  *  k  -  2.0  *  mu)  /  6.0; 
r3kd2  :=  (3.0  *  k)  /  2.0; 
sod3  ;=  so  /  3.0; 

Itimestep  :=  ntimesteps;  (last  time  step) 

dgamamax  :=  0.0;  (maximum  delta_gama  for  iteration) 

writelnCstart'); 

new(Fpiptr); 

new(sptr); 

new(Qptr); 

new(Tsvptr); 

new(nFpiptr); 

new(nsptr); 

new(nTsvptr); 

new(Tl  lavptr); 
new(T22avptr); 
new(T33avptr); 
new(T12avptr); 

writeCRestart  yy^n/ry  r"); 
ieadln(restait); 

write(TesiType  s/c/t  O; 
feadln(tesUype); 

wrile(TcstDim  f/r  :y, 
ieadln(testdim>; 

strainlimit  :=  'n'; 
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msdot  :=  1.0; 
if  testdim  =  'r'  then 
msdot  :=  -1.0; 


if  testtype  =  's'  then  begin 
open(Tfde,  TS.pr2'); 
open(pfile,  'p0s.pr2'); 
open(pOfile,  'pls.pr2'); 
open(p00file,  ’p2s.pr2’); 
end  {if} 

else  if  (tesUype  =  'c')  or  (testtype  =  ’f)  then  begin 
open(Tfde,  'TC.pr2'); 
open(pfile,  'p0c.pr2'); 
open(pOfile,  'plc.pr2'); 
open(p00file,  'p2c.pr2'); 
end; 

eulerread(totcrystals); 
genSc(Sc,  plllc); 
initialize; 

if  (restart  =  'nn')  or  (restart  =  'ry')  then  begin 
deltat  :=0.1;(2;} 

for  il  :=  2  to  ntimesteps  do  begin 
if  il  >  2  then 
deltat  :=  0.2; 
if  il  >  S  then 
deltat  :=  0.3; 
if  il  >  7  then 
deltat  ;=  0.4; 
if  il  >  10  then 
deltat  :=  0.6; 
if  il  >  15  then 
deltat  :=  1.0; 
if  il  >  20  then 
deltat  :=  2.0; 
if  il  >  30  then 
deltat  :=  4.0; 
if  il  >  50  then 
deltat  :=  7.0; 
if  il  >  80  then 
deltat  ;=  8.0; 
if  il  >  1(X)  then 
deltat  ;=  10.0; 
if  il  >  120  then 
deltat  :=  15.0; 

time[il]  :=  time(il  -  I]  +  (deltat / abs(straindot))  *  0.001; 
strain[il]  ;=  2.0  *  startstrain  -t-  msdot  *  straindot time(il]; 
(  writeln(time[il]);) 
end;(foril) 

writeln('Max.  Strain  =  ',  strain{ntimcstcps]); 


51 


end;  (if) 


for  il  :=  2  to  ntimesteps  do  begin 

10; 

if  (il  <=  10)  and  ((restart  =  'nn')  or  (restart  =  'ry'))  then  begin 
deltat  :=  time[il]  -  tiine[il  - 1]; 

strain[il]  2.0  *  startstiain  +  msdot  *  straindot  ♦  tiine[il]; 
end  (if) 
else  begin 

time[il]  :=  time[il  - 1]  +  deltat; 

strain[il]  :=  2.0  *  startstiain  +  msdot  *  straindot  *  timelil); 
end;  (else) 


ptsdgamamax  :=  dgatnamax; 
dgamamax  :=  0.0, 


si  :=  strain[il];  (straindot*  time[il];) 

if  testtype  =  's'  then  begin 
F  :=  133; 

F[1.2]  :=  si; 

end 

else  begin 

if  testtype  =  'c'  then  begin 
F  :=  null33; 

F[l.  1]  :=exp(sl/2.0); 

F(2.  2]  :=  F[l,  1); 

F(3, 3]  :=exp(-sl); 
end 

else  begin 
F  :=  null33; 

F[l.  l]:=exp(-sl/2.0); 

F[2,  2]  :=  F[l.  1]; 

F[3, 3]  :=exp(sl); 
end;  (else) 
end;  (else) 

mprint(F,  3); 

nTsvptr^  ;=  Tsvptr''; 
nsptr^  :=  sptr^; 


cdgamamax  :=  0.0, 

for  crystal  ;=  startcrystal  to  lastcrystal  do  begin 
writeln; 

writelnC**  Crystal  =  ',  crystal :  4, '  **'); 

for  alpha  :=  1  to  12  do  begin 
(write(’Sc',  alpha);) 

(mprint(Sc(alpha],  3);) 
aikbkj(Qptr^(crystal),  Sc  (alpha),  tcmp33); 

Sa(alpha)  ;=  null33; 
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for  i  :=  1  to  3  do 
for  j  :=  1  to  3  do 
for  1  :=  1  to  3  do 

Sa(alpha,  i.  j]  :=  Sa[alpha,  i.  j]  +  temp33[i,  1]  *  Qptr''[crystal,  j,  1]; 
syinmtopv(Sa[alpha],  Savec[^pha]); 

{ v6print(Savec[alpha]  ,6); ) 
end;  (for  alpha] 


writeln('Crystal=',  crystal :  4, ' :  time]',  il,  time[il] :  8  :  2,  ‘ :  strain[',  il,  strain[il]  :  8  :  3, ' :  dt=',  deltat 
8:4); 

aikbkj(F,  Fpipti^[crystal],  Feinter); 

{A  =  Transpose[Feinter]  J^einter) 

A  :=  null33; 
for  i  :=  1  to  3  do 
for  j  ;=  1  to  3  do 
for  1  ;=  1  to  3  do 

A[i,  j]  :=  A[i,  j]  +  Feinter[l,  i]  •  Fcintcr[l,  j); 

(  writeln('A  : ');] 

(  mprint(A,  3); ) 

si  :=  r3kin2inud6  *  tiace(A); 
s3  :=  si  -  r3kd2; 
for  i  :=  1  to  3  do 
for  j  ;=  1  to  3  do 

Tstr[i,  j]  :s  mu  *  A[i,  j]  +  s3  *  I33[i,  j]; 

{  mprint(Tstr,  3):} 

symmtov(Tstr,  Tstrvec); 

for  alpha  :=  1  to  12  do  begin 
Ba[^pha]  :=  null33; 
for  i  1  to  3  do 
for  j  :=  1  to  3  do 
for  I  :=  I  to  3  do 

Ba[alpha,  i,  j]  :=  Ba[alpha,  i,  j]  +  A[i,  1]  ♦  Sa[alpha,  1,  j]  +  Sa[alpha,  1,  i]  *  A[l,  j]; 

{  writeln('Ba[',alpha,T):) 

(  mprint(Ba[alpha],3);) 

si  ;=  trace(Ba[alpha])  •  r3km2mud6; 
for  i  :=  1  to  3  do 
for  j  :=  1  to  3  do 

Ca[alpha,  i.j]  ;=  mu  *  Ba[alpha,  i,j]  +  si  •  133[i.  j]: 

(  writeln('Ca[>lpha,']');) 

(  mprint(Ca[alpha],3);) 

symmtov(Ca[alpha],  Cavcc[alpha]): 
end;  (for  alpha) 
st  :=  nsptr^[crystal]; 

gamopdelt  :=  gamartoto  *  deltat; 
iterout  :s  false; 
iterate  1; 

(for  iterate  ;=  1  to  20  do  begin] 

repeat 
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tau_dgaina; 

{  vprint(tau.  12);) 

(  vprinKdeltagama,  12);} 

newdgamamax  ;=  absinax(deltagaina); 

{  writelnCndgamax,  dgamax  .ptsdgamax  = ',} 

{  newdgamamax:  10:4,' :  '.dgamamax:  10:4,' :  ',ptsdgamamax:  10:4); } 

if  (newdgamamax  >  dgamamax)  then 
dgam^nax  :=  newdgamamax; 
temp6  :s  nullv6; 
for  i  :=  1  to  6  do  begin 
for  alpha  :=  1  to  12  do 

temp6{i]  :=  temp6[i]  +  deltagama[alpha]  *  Cavec[alpha,  i); 

Tfii[i]  :=  nTsvptr^[cTystal][i]  -  Tstrvec[i]  +  temp6[i]; 
end;(for  i) 

writeln('C  = ',  crystal, ' :  Tstep  =  ',  il, ' :  iterate  =  ',  iterate); 

{  v6prini(Tfn,6);) 

Jacobian  :=  166; 
for  i  :=  1  to  6  do 
for  j  :=  1  to  6  do 
for  alpha  :=  1  to  12  do 

if  (abs(tau[alphal)  >  (so  /  100000.0))  then  begin 
si  :=  deltagama[alpha]  /  (m  *  taufalpha]); 

Jacobian[i,  j]  :=  Jacobian[i,  j]  +  si  *  Cavec[atpha,  i]  *  Savec[alpha,  j); 
end;  {if} 

(  writelnCJ ;  *);} 

(  m6print(Jacobian,6);} 

inverse6(Jacobian,  Uac,  6); 

(  writeln(U  : ');} 

(  m6print(Uac,Q;} 

temp6  :=  nullv6; 
for  i  :s  1  to  6  do  begin 
for  j  :=  1  to  6  do 

temp6[i}  :=  temp6[i}  +  Uac[i,  j}  •  Tfn[J}; 

Tsnewvec[i}  :=  nTsvptr^[crystal][i}  -  temp6[i}; 

Tvecerror[i}  ;=  Tsnewvec(i}  -  nTsvptr^[crystaI}[i}; 
end;(fori} 

Tout  :=  true; 

for  i  :=  1  to  6  do  begin 

if  ((abs(Tvecerror[i}))  >  so  /  100000.0)  then  begin 
Tout  :=  false; 

if  (abs(TveceiTor[i})  >  2.0  *  sod3)  then 
[  Tsnewvec[i}  :=  Tsnewvec[i}  +  nTsvptr^[crystal}[i})/2.0} 

(  else} 

Tsnewvec[i}  :=  nTsvpu^[crystal][i}  +  sod3  *  sign(TvecciTor(i}); 
end;  {if} 
end;  {fori} 

nTsvpir''(crystal}  :=  Tsncwvec; 

{ v6print(T  snewvec  ,6); } 

{  writeln(Tout :  ',Tout);} 

if  Tout  then  begin 
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tau.dgama; 
templ2  :=  nullvl2; 
for  alpha  :=  1  to  12  do  begia 
s3  :=  0.0; 

for  beta  :=  1  to  12  do  begin 
si  :=  1.0  -  nspti^[crystal,  beta]  /  ss; 
if  (si  >  0.0)  then  begin 
^  aa  •  ln(sl); 

si  :=  exp(s2);  (ho  embedded  in  qab) 

end 

else 

si  :=  0.0; 

templ2[alpha]  :=  templ2[alpha]  +  qab[alpha,  beta]  *  si  *  abs(deltagama[beta]); 
end;  (for  beta) 

snew[alpha]  :=  st[alpha]  +  templ2[alpha]; 
seiTor[alpha]  :=  snew[alpha]  -  nsptr^[crystal,  alpha]; 
end;  (for  alpha) 
nsptr^[ciystal]  :=  snew; 
sout  :=  true; 
si  :=  absmax(serror); 
if  si  >  so  /  100.0  then  begin 
sout  :=  false; 
end;  (if) 
write('s  O; 

vprint(nsptr^[crystal],  6); 
end;  {if  Tout] 

(  writelnCsout :  ’.sout);) 

iterate  ;=  iterate  +  1; 
if  (Tout  and  sout)  then 
iterout  :s  true; 

until  (iterout  or  (iterate  >  80)); 

{end;for  iterate) 

writeln('ndgamax,  dgamax  ,ptsdgamax  =  newdgamamax  :  10  :  4, ' :  dgamamax  :  10:4,':',  ptsdgamamax  :  10 
:4); 

if  (newdgamamax  >  cdgamamax)  then 
cdgamamax  :=  newdgamamax; 

{for  Mac) 

{  if  (il>10)  or  (restartss'y')  then) 

{  if  (iterate>=60)  or  (dgamamax>  100000000.0)  then  begin] 

{  writelnC*’*""*  deltat  is  large  -  unstable  -  reduce!!  **♦♦■);) 

{  deltat  :=  0.7S'*deltat;} 

{  goto  10:) 

{  end;) 

(if) 


temp233  :=  133; 
for  alpha  :=  1  to  12  do  begin 
smuim(dellagama{alpha],  Sa{alpha].  tempi  33); 
for  i  :=  1  to  3  do 
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for  j  :=  1  to  3  do 

temp233(i,  j]  :=  temp233[i.  j]  -  templ33[i,  j]; 
end;(fora^ha) 

aikbkj(Fpiptr'(crystal],  temp233,  tempi  33); 
si  :=  det33(templ33); 
s2  :=  exp(ln(sl)  /  3.0); 
si  :=  l/s2; 

smulm(sl,  templ33,  Fpinew); 
writelnCFpinew :  ■);) 

mprint(Fpiiiew3): ) 

nFpipir^[crystal]  :=  Fpinew; 

aikbkj(F,  Fpinew,  Fe); 
writelnCFe :  *);) 
mprint(Fc,3);) 

vtosymm(nTsvptr*[crystal],  Ts); 
aikbkj(Fe,  Ts,  temp33); 
si  ;=  det33(Fe); 

Tcauchy  :=  null33; 
for  i  :=  1  to  3  do 
for  j  :=  1  to  3  do  begin 
for  1  :=  1  to  3  do 

Tcauchy[i,  j]  :=  icauchyii,  j]  +  temp33[i,  1]  *  Fe(j,  1); 

Tcauchy[i,  j]  ;=  Tcauchy [i,j]  /  si; 
end;{forj} 
writeln(Tcauchy :  O: 
mprint(Tcauchy,  3); 

if  (crystal  s  staitcrystal)  then  begin 
T12avptr^(il]  :=  Tcauchy(l,  2]; 

Tllavptr^[il]  ;=Tcauchy[l,  1]; 

T22avptr''[il]  :*  Tcauchy[2, 2); 

T33avptr^[il]  :=  Tcauchy[3, 3]; 
end  (if) 
else  begin 

itemp  :=  crystal  -  staitcrystal  +  1; 

T12avptr'[il]  :=  T12avptr''[il]  +  (Tcauchy [1,2]  -T12avptr^[il])/ itemp 
Tllavptr''[il]  :=  Tllavptr''[il]  +  (Tcauchy[l,  1]  -  Tllavptr''{il])/  itemp 
T22avptr^[il]  :=T22avptr''[il]  +  (Tcauchy[2, 2]  -T22avptr''[il])/ itemp 
T33avptr''[il]  :=  T33avptr''[il]  +  (Tcauchy [3, 3]  -  T33avptr^[il])  /  itemp 
end;  [else) 

if  testdim  =  T  then 
if  (strain[il]  >  maxstiain)  then 
strainlimit  :=  'y'; 
if  testdim  =  'r'  then 
if  (strain[il]  <  minstrain)  then 
strainlimit  :=  'y'; 

if  (il  =  ntimesteps)  or  (strainlimit  =  'y')  then  begin 
RUdecomp(Fe,  R,  InU); 
aikbkj(R,  Qptr''[crystal],  temp33); 


writelnCpfile, "); 
for  i  :=  1  to  8  do  begio 
indoiv(temp33,  pi  1  lc[i],  newpoie); 
if  (newpole[3]  >=  0.0)  then  begin 
eapolexy(newpole,  x,  y); 
write(pWe,  x  :  10  :  4,  y  :  10  :  4); 
end;  (if) 
end;  (for  i) 

writeln(pOfilc, "); 
for  i  :=  1  to  12  do  begin 
mdotv(temp33,  pi  10c[i],  newpoie); 
if  (newpole[31  >=  0.0)  then  begin 
eapolexy(newpole,  x,  y); 
write(p0file,  x  :  10  :  4,  y  :  10  : 4); 
end;  (if) 
end;(fori} 

writeln(p00file, "); 
for  i  :=  1  to  6  do  begin 
mdotv(temp33,  pl00c[i],  newpoie); 
if  (newpole[3]  >=  0.0)  then  begin 
eapolexy(newpole,  x,  y); 
write(p00file.  x  :  10  :  4,  y  :  10  :  4); 
end;  (if) 
end;  (fori) 

end;  (if) 

end;  (for  crystal) 


writeln('dt=  deltat :  8  :  3, '  timeC,  il,  ']=  ume[il]  :  8  :  2, '  strain=',  strain[il]  :  8  :  4) 
Rt  :=  edgamamax  /  sdgama; 

(auto  time  increment  selection) 
if  dgamamax  <  1000000.0  then  begin 
if  Rt  <=  0.8  then 
deltat  :=  1.2S  *  deltat 
else  if  (Rt  >  0.8)  and  (Rt  <  1.25)  then 
deltat  :=  deltat  /  Rt 
else  begin 

writelnC****  deltat  is  large  -  reduce!  ****’); 
deltat  :=  0.75  •  deltat; 
end;  (else) 
end  (if) 
else  begin 

writelnC****  deltat  is  large  -  unstable  -  reduce!!  *****); 
deltat  :=  0.75  *  deltat; 
end;  (else) 

writelnCcdgamamax,  Rt edgamamax  :  10  :  4,  Rt :  10  ;  4); 
writeln('newdeltat  = deltat :  10  : 4); 
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Fpipcr^  :=  nFp4)tr^; 
sptr^  :=  nsptr^; 

Tsvpii^  :=  nTsvptt^; 

writelnCTfile,  strain[il] :  12:4,  T12avptr^[il] :  12  :  4,  T33avptr^[il] :  12:  4,  T22avptr''[il] :  12  :  4,  TI  lavptr''[il] 
12 :  4.  -<T33avpti^[il]  -  (Tllavpti^[il]  +  T22avptr^[il])/2.0) :  12  :  4); 

if  strainlimit  =  'y'  then  begin 
Itimestep  :=  il; 
goto  20; 
end;  (if) 

end;(far  il) 

20: 

writeln(il  :  5,  ctystal :  5); 

close(Tfile); 

cIose(pfilc); 

close(pOfile); 

close(p00flle); 

saversfiles; 

end. 
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